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ABSTRACT.  Machine  tool  chatter  has  been  characterized  as  isolated  periodic  solutions  or  limit  cycles  of 
delay  differential  equations.  Determining  the  amplitude  and  frequency  of  the  limit  cycle  is  sometimes  crucial 
to  understanding  and  controlling  the  stability  of  machining  operations.  In  Gilsinn  [9]  a result  was  proven 
that  says  that,  given  an  approximate  periodic  solution  and  frequency  of  an  autonomous  delay  differential 
equation  that  satisfies  a certain  non-criticality  condition,  there  is  an  exact  periodic  solution  and  frequency  in 
a computable  neighborhood  of  the  approximate  solution  and  frequency.  The  proof  required  the  estimation  of 
a number  of  parameters  and  the  verification  of  three  inequalities.  In  this  paper  the  details  of  the  algorithms 
will  be  given  for  estimating  the  parameters  required  to  verify  the  inequalities  and  to  compute  the  final 
approximation  errors.  An  application  will  be  given  to  a Van  der  Pol  oscillator  with  delay  in  the  nonlinear 
terms.  A MATLAB  m-file  implementing  the  algorithms  discussed  in  the  paper  is  given  in  the  appendix. 
AMS  (MOS)  Subject  Classification.  34K11,  34K13,  34K28. 

1.  Introduction 

Machine  tool  dynamics  has  been  modeled  using  delay  differential  equations  for  a number 
of  years  as  is  clear  from  the  vast  literature  associated  with  it.  For  a detailed  review  of 
machining  dynamics  see  Tlusty  [32].  For  a discussion  of  dynamics  in  milling  operations  see 
Balchandran  [1]  and  Zhao  and  Balachandran  [34],  For  drilling  operations  see  Stone  and 
Askari  [29]  and  Stone  and  Campbell  [30].  For  an  analysis  of  chatter  occurring  in  turning 
operations  see  Hanna  and  Tobias  [16],  Marsh  et  al.  [21],  and  Nayfeh  et  al.  [22].  Machine 
tool  chatter  is  undesirable  self-exited  periodic  oscillations  during  machining  operations.  It 
has  been  identified  as  a Hopf  bifurcation  of  limit  cycles  from  steady  state  solutions.  For  a 
way  of  estimating  the  critical  Hopf  bifurcation  parameters  that  lead  to  machine  tool  chatter 
see  Gilsinn  [8]. 

In  studying  the  effects  of  chatter  it  is  sometimes  desirable  to  compute  the  amplitude  and 
frequency  of  the  limit  cycle  generating  the  chatter.  This  entails  solving  the  delay  differential 
equations  that  model  the  machine  tool  dynamics.  There  is  a large  literature  on  numerically 
solving  delay  differential  equations.  Some  representative  methods  are  described  in  Banks 
and  Kappel  [2],  Engelborghs  et  al.  [5],  Kemper  [20],  Paul  [23],  Shampine  and  Thompson 
[25],  and  Wille  and  Baker  [35].  Although  these  methods  generate  solution  vectors  that  can 
be  studied  by  harmonic  and  power  spectral  methods  to  estimate  the  frequency  of  periodic 
cycles,  they  do  not  directly  generate  a representative  model  of  a limit  cycle  such  as  a Fourier 
series  representation. 

It  is  also  desirable  to  know  whether  a representation  of  an  approximate  limit  cycle  is 
close  to  a true  limit  cycle.  In  other  words  we  wish  to  answer  the  question  as  to  whether  the 
approximate  solution  represents  sufficiently  well  a true  solution.  This  is  answered  with  a set 
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of  test  criteria  by  Gilsinn  [9],  who  showed  that,  given  a representative  approximate  solution 
and  frequency  for  a periodic  solution  to  the  autonomous  delay  differential  equation 

(1)  x = X.(x(t),x(t  - h)), 

where  x,X  £ Rn,  h > 0,  X sufficiently  differentiable,  there  are  conditions,  depending  on 
a number  of  parameters,  for  which  (1)  has  a unique  exact  periodic  solution  and  frequency 
in  a computable  neighborhood  of  the  approximate  solution  and  frequency.  This  result  was 
first  established  in  a very  general  manner  for  functional  differential  equations  by  Stokes  [28] 
who  extended  an  earlier  result  for  ordinary  differrential  equations  in  Stokes  [27].  A crucial 
aspect  in  applying  the  result  involves  verify  a certain  ” non-criticality”  condition.  However, 
no  computable  algorithms  were  given  in  the  case  of  functional  differential  equations  to 
estimate  the  various  parameters  or  verify  the  ” non-criticality”  condition.  Only  recently 
have  algorithms  been  developed  to  computationally  verify  these  conditions  in  the  fixed 
delay  case.  A preliminary  announcement  of  algorithms  for  computing  these  parameters  and 
verifying  the  ” non-criticality”  condition  was  given  by  Gilsinn  [7].  In  this  paper  we  include  a 
more  detailed  discussion  of  the  algorithms  and  apply  them  to  a Van  der  Pol  equation  with 
delay  in  its  nonlinear  terms. 

The  result  of  Stokes  [28]  for  functional  differential  equations  depends  on  verifying  certain 
conditions  that  require  computing  various  parameters,  hi  order  to  apply  Stokes’  result 
a proof  in  the  case  of  equation  (1)  will  be  given  here  since  certain  inequalities  that  are 
developed  within  the  proof  are  necessary  for  proving  the  fixed  point  contraction  mapping 
conditions  and  rely  on  the  specific  form  of  (1). 

The  notation  used  in  the  paper  is  described  in  Section  2.  The  non-criticality  condition  is 
defined  in  Section  3.  In  Section  4 we  construct  an  exact  frequency  and  27r-periodic  solution 
of  (1)  as  a -perturbation  problem.  In  Section  5 we  define  a map  that  is  used  to  prove,  by  a 
contraction  argument,  the  existence  of  an  exact  frequency  and  27r-periodic  solution  of  (1). 
The  main  contraction  theorem  is  proven  in  Section  6.  A Galerkin  algorithm  to  compute  a 27r- 
periodic  solution  to  a nonlinear  autonomous  delay  differential  equation  is  given  in  Section  7. 
The  Floquet  Theory  for  delay  differential  equations  is  discussed  in  Section  8.  An  algorithm 
for  computing  the  characteristic  multipliers  of  the  variational  equation  of  (1)  with  respect  to 
the  approximate  27r-periodic  solution,  is  outlined  in  Section  9.  An  algorithm  to  determe  the 
solution  to  the  formal  adjoint  equation  with  respect  to  the  variational  equation  of  (1)  with 
respect  to  the  approximate  27r-periodic  solution,  is  outlined  in  Section  10.  An  algorithm  for 
estimating  a critical  parameter  is  given  in  Section  1 1 . An  application  of  these  algorithms  to 
the  Van  der  Pol  equation  with  delay  is  given  in  Section  12.  Conclusions  are  given  in  Section 
13  and  a disclaimer  is  given  in  Section  14.  The  derivation  of  the  differentiation  matrix  (86) 
is  given  in  Appendix  1,  Section  15.  Certain  bounds  and  Lipschitz  conditions  used  in  the 
fixed  point  theorem  are  proven  in  Appendices  2 and  3 (Section  16  and  Section  17).  The 
MATLAB  code  and  associated  support  functions  implementing  the  algorithms  are  given  in 
Appendix  4,  Section  18. 


2.  Notation 

Let  CL  denote  the  space  of  continuous  functions  from  [— u,  0]  to  C"  with  norm  in  C ^ 
given  by  \(f)\  — max|<?!>(s)|  for  — u>  < s < 0,  where 

(2)  \Hs)\=  (j2\Us)\2 

C L is  a Banach  space  with  respect  to  this  norm.  Let  V be  the  space  of  continuous  27r-periodic 
functions  with  sup  norm,  |-|  on  (—  oo,  oo).  Let  V\  C V be  the  subspace  of  continuously  differ- 
entiable 27r-periodic  functions  with  the  sup  norm.  Let  X(x,  y)  be  continuously  differentiable 
in  some  domain  Qn  C Cn  x Cn  with  bounded  derivatives  where 

Xi(x,y)\<B, 


(3) 
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for  i = 1.2,  (x,y)  € f2n.  The  subscripts  of  X indicate  derivatives  with  respect  to  the  first 
and  second  variables  of  X respectively.  We  further  assume  that  the  first  partial  derivatives 
satisfy  Lipschitz  conditions  given  by 

(4)  \Xi(xi,yi)  - Xl(x2,y2)\  < IC(\xi  - x2\  + \yi  - y2\), 
for  {xi,yi),  (x2,y2)  € &n- 

In  order  to  simphfy  the  notation  for  (1)  we  will  first  normalize  the  delay  h to  unity  by 
setting  s — t/h.  Then.  (1)  becomes 

(5)  j(s)  = /iX(y(s),y(s-l)) 

as 

where  y(s)  = x(sh).  Therefore  we  will  assume  h = 1 in  (1).  We  will  also  make  one  further 
transformation.  Since  the  period  T = 2t /uj  of  a periodic  solution  for  (1)  is  unknown  we  can 
normalize  the  period  to  [0.  27t]  by  introducing  the  substitution  of  t/uj  for  t and  rewriting 
(1).  with  h = 1,  in  the  form 

(6)  uj  x = X(x(t),  x(t  — uj)). 

For  ipi,ip2  € V we  denote  the  total  derivative  of  X(x,y)  by 

(7)  dX{x,y,%l)i,ip2)  = Xi{x,y)ipi  + X2(x,y)jjj2. 

Let  A(t),  B(t)  be  continuous  27r-periodic  matrices.  Then  a characteristic  multiplier  is 
defined  as  follows. 


Definition  2.1.  p is  a characteristic  multiplier  of 

(8)  y = A(t)y(t)  + B(t)y(t  — u) 

if  there  is  a nontrivial  solution  y(t)  of  (8)  such  that  y(t  + 2ir)  = p y(t).  Note  that  if  p = 1 
then  y(t)  is  2n-periodic. 


To  simplify  some  of  the  notation  we  will  suppress  the  t and  write,  for  example,  x = 
x(t),  x u = x(t  — uj),  but  in  other  cases  we  will  maintain  the  t,  especially  when  describing 
computational  steps.  We  will  also  at  times  use  the  notation 

n 1/2 


(9) 


F|2  = 


/*  27T 

Jo 


|x(f)|"  dt 


3.  Non-criticality  Condition 


Galerkin  and  harmonic  balance  methods  can  be  used  to  develop  27r-periodic  approximate 
solutions  for  (6).  A fast  discrete  Fourier  series  algorithm  for  computing  an  approximate 
series  solution  and  frequency,  (lj,x),  has  been  given  by  Gilsinn  [8].  See  Section  7 below  for 
a brief  discussion  of  a Galerkin  method  for  approximating  a solution.  At  this  point,  then, 
we  assume  that  we  have  developed  a 27r-periodic  approximate  solution  and  frequency,  (uj,  x) 
for  (6),  where  x is  27r-periodic  and 

(10)  ujx  = X(x,  xcu)  + k 
where  k(t)  is  a 27r-periodic  residual  bounded  by 

(11)  |fc|  < r. 

The  required  size  of  the  residual  error,  r,  will  become  clear  based  upon  estimates  later  in  this 
paper.  These  estimates  will  indicate  in  particular  situations  how  accurately  an  approximate 
solution  and  frequency  would  need  to  be  computed. 

The  variational  equation  with  respect  to  the  approximate  solution  and  frequency  is  given 
by 

(12)  Cjz  = dX  (x,xoj\z,  Z&) . 

Let  A = X\  (x,x;f),  B = X2  (x,xcj).  The  formal  adjoint  of  (12)  is  given  in  row  form  by 

(13)  ujv  = — vA  — V-cjB. 
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The  next  lemma,  proven  in  Halanay  [14],  relates  the  number  of  independent  27r-periodic 
solutions  of  (12)  to  those  of  (13). 

Lemma  3.1.  System  (12)  and  (13)  have  the  same  finite  number  of  independent  2n -periodic 
solutions. 


We  will  not  give  the  proof  of  the  next  lemma,  since  it  is  also  proven  in  Halanay  [14].  The 
result,  however,  will  be  critical  to  the  main  approximation  theorem. 

Lemma  3.2.  The  nonhomogeneous  system 

(14)  lox  — Ax  + Bxcj  + / 
has  a unique  2n -periodic  solution  if  and  only  if 

p27T 

(15)  / vlfdt  = 0 

Jo 

for  all  independent  solutions  Vo  of  period  2ir  of  (13).  Furthermore  there  exists  an  M > 0, 
independent  of  f,  such  that 

(16)  \x\<M\f\. 

We  will  give  the  proof  of  the  next  lemma.  Although  it  is  stated  in  Hale  [15]  and  in  Halanay 
[14],  the  proof  is  not  generally  available.  The  result,  however,  motivates  the  definition  of  a 
non-critical  approximate  solution. 


Lemma  3.3.  Let  p = 1 be  a simple  (i.e.  multiplicity  one)  characteristic  multiplier  of  (12) 
Let  p be  a non-trivial  solution  of  (12)  associated  with  p.  Define 

(17)  =p+  Bpz, 
then 

/•27T 

(18)  / voTJ(p, u)dt  0 

Jo 

for  all  independent  solutions  v of  the  adjoint  (13). 

Proof.  Let  y(t)  be  any  27r-periodic  solution  of  (12)  and  write 

(19)  z = y + tp 
Then,  substituting  (19)  into  (12),  write 

Coz  = Az  + Bzjj  + Si  ^ p + Bpx) 

(20)  = Az  + Bz^j  + ujJ  ( p . lu) 

We  will  suppose 

r2n 

(21)  / vT  J{p,La)dt  = 0 

Jo 

and  show  a contradiction.  Since  p is  a simple  characteristic  multiplier  then  p is  the  only 
non-trivial  27r-periodic  solution  associated  with  p and  so  there  is  only  one  solution  of  (13) 
associated  with  1/p.  Lemma  3.2  and  (21)  imply  that  there  is  a unique  27r-periodic  solution 
2 of  (20).  With  2 and  p both  27r-periodic  let  w — z — tp.  w cannot  be  a multiple  of  p for 
there  would  be  a to  such  that  top  = 2 — tp  or  2 = (to  + t)  p.  Since  z,p  are  27r-periodic  we 
have  2 — (to  + t + 2tt)  p.  But  then  we  would  have  to  — to  + 2i r or  2ir  = 0,  a contradiction. 
QED 

Lemma  3.2  will  imply,  in  the  case  of  a non-critical  27r-periodic  approximate  solution  of 
(6).  that  there  is  only  one  Vq  in  (22). 

We  can  now  give  the  definition  of  a non-critical  approximate  solution  of  (6). 
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Definition  3.4.  The  pair  (Co,x),  where  x is  at  least  twice  continuously  differentiable,  is  said 
to  be  non-critical  with  respect  to  (6)  if  (1)  the  variational  equation  about  the  approximate 
solution  x,  given  by  (12),  has  a simple  characteristic  multiplier  po  with  all  of  the  other 
characteristic  multipliers  not  equal  to  one.  (2)  If  Vq,  l^ob  = 1,  is  the  solution  of  (13) 
corresponding  to  pQ,  i.e.  with  multiplier  1/po,  then 


(22) 

where 

(23) 


/ 


vq 1 J(x,Co)dt  ^ 0, 


J(x,Co)  = x + Bx& 


4.  A Perturbation  Problem 


In  this  paper  we  will  look  for  an  exact  27r-periodic  solution,  x,  and  an  exact  frequency, 
uj,  for  (6)  as  a perturbation  of  the  27r-periodic  approximate  solution  , x , and  approximate 
frequency,  Co,  of  (6).  In  particular,  let 


uo  = Co  + (3 


(24) 


X = X -\ 2 

CO 


Then,  substituting  (24)  into  (6)  and  using  (10),  we  can  write  the  equation  for  2 and  (3  as 

(25)  luz  = dX  ( x , x^;  z , z&)  + R{z,(3)  — (3J  (x,  Cj  j — k 
where 

(26)  R(z,(3)=  X [x + —z,xw + —zj\  — X (x,xc, 

\ uj  u 


- dX  (x,  xa>',z,  zcf)  + (3Bxcj. 


and  J is  given  by  (23). 

In  the  next  lemma  we  establish  bounds  and  Lipschitz  conditions  for  R{z,(3). 


Lemma  4.1.  There  exist  functions  TZo(z,(3)  > 0.  lZi{z,  (3,z,  (3)  >0,  i = 1,2,  such  that 
IZq  — > 0 as  (z,  /?)—>•  0 and  1Z,  — > 0 os  (z,  (3 , z,  f3)  — > 0 and 


(27)  \R(z,(3)\  < n0(z,(3) 

\r(z,0)-R(zJ)\  < 'Jl1(z,p,zj)\z-z\  + n2(z,p,zj)\(3-p\ 

Proof:  Appendix  2. 

Since  we  will  be  considering  \(3\  small,  we  will  begin  by  restricting  f3,  which  could  be 
negative,  so  that 


(28) 


Co  + (3  > 


UJ 

9 ' 


We  can  select  \(3\  < Co/2. 

As  a first  step  to  establishing  the  existence  of  a 2/r-periodic  solution  of  (25)  we  first  study 
the  existence  of  a 27r-periodic  solution  of 


(29)  Coz  = dX  (x,  xjj ; z,  z^)  + g — (3J  ^x,  u?J  — k 

where  g G V.  For  this  we  have  the  following  lemma 


Lemma  4.2.  If  (Co,x)  are  non-critical  with  respect  to  (6),  then  (a)  there  exists  a unique  (3 
such  that 

(30)  g — (3J  ^x,u;j  — k Cl  vo 

where  vq  is  the  solution  of  (13)  corresponding  to  the  characteristic  multiplier  po  of  (12),  and 
(b)  there  exists  a unique  2n -periodic  solution  of  (29)  that  satisfies 
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for  some  M > 0. 


Proof:  Take 


(32) 
where 

(33) 


and  apply  Lemma  3.2. 

We  can  now  establish  bounds  on  0,  z and  i.  For  notation,  designate  the  unique  0 and 
2 in  Lemma  4.2  by  0(g)  and  z(g)  respectively,  and  z by  z(g). 


Lemma  4.3.  There  exist  three  constants,  designated  by  A,,  i — 0,  1,  2,  such  that 

\0(g)\  < AoGsl  +0 

(34)  \z(g)\  < Ai(|-p|  + r) 

\z(g)\  < A2(M  + r) 


Proof:  From 

I5I2  < V2rrr  |(/| 

(35)  \k\2  < V2^  |&| 
and  the  Cauchy-Schwarz  inequality  applied  to  (32) 

(36)  \0(g)\  < |a|  |v0T|2  | g - k\2  < V2n\a\(\g\  + r) 

from  the  bound  \k\  < r. 

From  (31) 


z(g) | < M [|p|  + \k\  + \0(g)\  | J (F,ch)  |]  , 


(37) 

< AI  \ 1 + V2n  ct| 

J ^x,  (2^ 

] (lsl  + 0 

From  (29) 

H|i(s)|  < 

1 dX  (x,  x^;z(g),  z(g)^)\  + \g 

— 0J  ^x,<hj  — A- 

(38) 

< 

[B(\z(g)\  + k(5)<o|)] 

(39) 

< 

(2MB)  [l  + y/2ir\a\ 

J ^x,u;^ 

] (Isl  + 0 

Therefore,  from  (36),  (37),  (38), 

Ao  = \/27r|a|, 

(40)  Ai  = M j^l  + v^|o|  | J |J 

A2  = -4t7'1  + 2 MB) 


5.  A Map  and  its  Properties 


In  the  main  approximation  theorem  we  will  show  that  the  solution  of  the  perturbation 
problem  (25)  is  the  fixed  point  of  a particular  contraction  map.  In  this  section  we  will  define 
the  map  and  establish  some  properties. 

We  begin  by  defining  a subset  of  V,  designated  by  Ms,  as 

Ms  = {g  e V : \g\  < <5}, 


(41) 
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where  <5  > 0.  Following  Stokes  [28]  we  will  define  a map  S : Ms  — >■  V in  terms  of  two 
mappings 

L : Ms^RxV  i, 

(42)  T : RxV  i~>V. 

To  define  L,  let  g G Ms,  then  Lemma  4.2  assures  us  of  the  existence  of  a unique  (3(g) 
satisfying  (30)  and  a unique  solution  z(g ) satisfying  (29).  Thus,  define  L : Ms  — >•  R x V\  by 

(43)  L(g)  = (0(g),  z(g)). 

Now  define  T : R x V\  — > V by 

(44)  T(0,  z)  = R(z,  0). 

Finally,  define  S : Ms  — > P by 

(45)  S(g)  = T(L(g))  = R(z(g),0(g)). 

Lemma  5.1.  For  g € A/5,  5 G Ms  there  exist  two  functions  Ei(5),  E2(S)  and  two  positive 
constants  F\,  F2  so  that 


(46) 

|S(S)|  < 

Ex(5), 

\S(g)-  5(9)|  < 

E2(8)  \g-g\, 

where 

(47) 

E08)  < 

F182, 

E2(5)  < 

f2s. 

Proof: 

From  (45)  and  (27)  we  have 

(48) 

\S(g)\  < TZ0  (z(g),0(g)) 

\S(g)~S(g)\  < nl(z(g),0(g), 

z (9),  0(g))  \z(g)  - 2 

(5)1 

+R 2 (z(g),  0(g),  z (g) , 0(g)) 

1 0(g) 

-0(g)  1 

By  Cauchy- Schwarz,  the  fact  that  \vqT | = 1,  and  (40) 

/ 

2tt 

\0(g)-0(g)\  < M/ 

J 0 

KT  (5  - 9) 

dt, 

r27T 

1/2 

(49) 

< |a| 

/ \g-g\~dt 
0 

’ 

< ^0 1 g 

-g\  ■ 

From  Lemma  4.2  and  the  definition  of  0(g), 

0(g)  we  have 

f27T  rr,  / / . 

\ \ 

J v0T  (g  - 0(g)J  (x,u 

j)  — k J dt  — 0 

(50) 

J voT  [g  - 0(g)  J (x,i 

— k^j  dt  = 0 

Then,  by 

subtracting. 

r2n 

(51) 

J vo  [(g  - g)  - (0(g)  - 0 (g))  J 

dt  = 

0. 

Lemma  4.2  also  shows  that  there  exists  a unique  z such  that  there  exists  a z such  that 
(52)  d)z  = dX  (x,x&-,z,zcj)  + [(3  -9)-  (0(g)  ~ 0 (9))  J . 

But  from  (29),  z(g)  — z (g)  also  satisfies  (52),  so  that  z = z(g)  — z (g)  and  from  (31) 
\z{g)-z(9) I < M | (5  -g)  - {0(g)  - 0(g))  J (x,u)  | , 

< M\g-g\- 


(53) 
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Then  (46)  follows  from  (48)  through  (53).  QED 

The  specific  forms  for  E\(8)  and  FA^)  are  given  in  Appendix  A. 2 as  equations  (201)  and 
(204),  respectively,  as  well  as  the  selection  of  F\  and  F2  as  (202)  and  (205)  respectively. 
As  functions  of  the  other  parameters  Ei(8)  and  depend  linearly  on  K,  and  B.  but 

non-linearly  on  Aq,  Ai,  and  A2  and  thus  non-linearly  on  M. 


6.  Main  Approximation  Theorem 


In  the  main  theorem  the  constants  F\,  F2  are  those  from  Lemma  4.2. 


Theorem  6.1.  If  (a)  (Co,x)  is  non-critical  with  respect  to  (6)  in  the  sense  of  Definition  3-4,. 
(b)  8 is  selected  so  that 

(54)  8 < minjl/Ei,  I/2F2,  lu/A\q} 

and  (c)  r < 8,  then  there  exists  and  exact  frequency,  u >* , and  solution,  x* , of  (6)  such  that 

\x*  — x\  < 4Ai<5, 

(55)  \u*  — o)|  < 2A  08, 

where  Aq,Ai  are  defined  in  (40)  and  8 is  defined  in  (41)- 


Proof:  Let  (3  and  2 be  defined  as  in  (24).  By  substituting  (24)  into  (6)  we  have  (25). 
Associated  with  (25)  we  consider  (29).  We  then  define  the  set  Afs  in  (41)  and  consider  the 
map  S : A f$  — > V defined  in  (45).  From  (200),  (201),  and  (202)  we  have  |5(g)|  < F\82 
for  g € Afs-  Furthermore,  from  (203),  (204),  and  (205)  we  have,  for  g,  g E Afs,  that 
| S(g)  — 5(p)|  < F28  | g — g |.  Now,  if  we  select  8 as  in  (54)  then  F\82  < 8 and  F28  < 1/2,  S 
maps  Afs  to  itself  and  is  a contraction.  The  last  inequality  that  8 satisfies  in  (54)  assures 
that  0(g)  satisfies  (28)  by  way  of  Lemma  4.3,  provided  r satisfies  r < 8.  Therefore,  S has 
a fixed  point  g*  £ Afs ■ This  implies  that  there  exists  a unique  (0*,z*),  z*  is  27r-periodic, 
satisfying  (25).  Then,  from  (24),  there  exists  a unique  (lo*,x*),  x*  is  27r-periodic,  satisfying 
(6).  From  (24),  with  r < 8, 


|u 1*  - Co\  < \0*\  < A0  (|<7*|  + r)  < 2 A0<5, 

2*1  < 2Ax  (\g*\+r)  < 4Ai<5. 

QED 

We  need  to  introduce  a note  here  on  the  relationship  between  r and  8.  I11  practice  the 
process  of  determining  them  is  iterative.  We  start,  by  determining  an  approximate  solution 
and  the  residual  r.  We  then  compute  all  of  the  parameters  that  involve  the  approximate 
solution  and  compute  8 from  (54).  We  then  compare  r against  8.  If  r < 8 we  are  finished, 
otherwise  we  have  to  return  and  recompute  another  approximate  solution  with  possibly 
smaller  residual  r and  iterate  the  process.  The  author  is  not  familiar  with  any  result  that 
guarantees  that  at  some  point  r < 8,  although  he  suspects  that  this  will  eventually  happen 
in  most  practical  problems. 


(56) 


\x*  — x I < 


7.  Approximating  a Solution  and  Frequency 

An  approximate  solution  and  frequency  for  (6)  can  be  developed  by  assuming  a finite 
trigonometric  polynomial  of  the  form 

m 

(57)  xm  = 02  cos  t + ^ [a2n  cos  nt  + 02n-i  sinnt] 

n= 2 

where  the  sinf  term  has  been  dropped  so  that  we  can  estimate  a\  = Co,  the  frequency. 
Note  that  we  have  centered  the  approximate  solution  about  the  origin,  since  we  assumed 
AQO.  0)  = 0.  If  we  set  a = (01,  02,  • • • , 02 m),  and 

Emit,  a)  = a\Xm(t)  - x ( Xm(t ),  Xm(t  ~ Oi)) 


(58) 
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then  for  a sufficiently  fine  mesh,  specified  by  {tt  : i = 1.  2,  • • • , 2 TV},  in  [0.  27t],  where 

2 i - 1 

(59)  t,  = —n, 

the  determining  equations  for  a can  be  written  as  (see  Urabe  and  Reiter  [33]) 


(60) 


*i(a)  = 


F2(  a)  = 


^2n-l(a)  = 


^2n(a)  = 


iV 

J_ 

Iv 

1 

Iv 

1 

N 


2N 

y Em  (it,  a)  sin tj  = 0 

i=i 

2N 

{ti,  a)  cos  ti  — 0 

i= 1 
2N 

-Em  (ti,  a)  sin  nti  = 0 

i=l 

2N 

Em  (tj,  a)  cos  nti  — 0 

i=l 


for  n = 2,  ■ ■ ■ , m. 

These  equations  give  2m  equations  in  2m  unknowns.  Standard  numerical  solvers,  using, 
for  example,  Newton’s  method,  for  nonlinear  equations  can  be  used  to  solve  for  a.  The 
number  of  harmonics,  m,  and  the  quadrature  index,  N.  can  be  selected  independently. 


8.  Floquet  Theory  for  DDEs 

The  analysis  of  the  stability  of  an  approximate  periodic  solution  for  (1)  usually  involves 
the  following  considerations.  If  x(t),  x £ Rn  is  an  approximate  periodic  solution  of  (1)  of 
period  2tt.  and  Q an  approximate  frequency,  then  the  linear  variational  equation  about  x(t) 
can  be  written 

(61)  z(t)  = A(t)z(t)  + B(t  )z(t  — 0), 

where  A(t)  and  B(t ) were  defined  previously  and  are  periodic,  with  period  2tt. 

We  now  define  the  period  map  U : Cb  — » Co  with  respect  to  (61)  by 

(62)  (U 4>)  (s)  = z(s  + 2tt) 

where  z(s)  is  a solution  of  (61)  satisfying  z(s)  = </>(s)  for  s £ [ — tD,  0] . In  this  paper  we 
assume  u < 2tt.  U is  then  a compact  operator  on  Cq,  whose  spectrum  is  at  most  countable 
with  0 as  the  only  possible  limit  point  (Halanay  [14]). 

A Floquet  theory  for  (61)  has  been  developed  by  Stokes  [26]  . In  particular,  if  a(U) 
represents  the  spectrum  of  U,  then  for  each  A £ a(U),  U<j)  = \<j>.  That  is,  the  spectrum 
consists  of  eigenvalues.  Furthermore,  the  space  Cq  can  be  decomposed  as  the  direct  sum  of 
two  invariant  subspaces 

(63)  Cb  = £(A)©A'(A) 

E(  A)  is  finite  dimensional  and  composed  of  the  eigenvectors  with  respect  to  A.  Furthermore, 
a {U\k)  = & (U)  — {A}.  If  {ipi},i  — 1,  • • ■ , d is  a basis  for  E( A)  and  we  let  T be  the  matrix 
with  columns  tpj  for  j = 1,  ■ ■ • , d,  then  there  is  a matrix  G(A)  such  that 

(64)  ET  = TG(A) 

Thus  we  can  think  of  Cq  as  being  a countable  direct  sum  of  the  invariant  subspaces  E(\,) 
plus  a possible  remainder  subspace,  B.  That  is 

(65)  Cq  = E(A1)©E(A2)©---0i? 

where  R is  a ’’remainder”  set  in  which  any  solution  of  (61)  with  initial  condition  in  R decays 
faster  than  any  exponential. 

For  each  of  the  E(Xl)  there  is  a basis  set  T,,  and  a matrix  G{ \j).  If  we  define  an  at  most 
countable  basis  set  {'F,},  i = 1,  2,  • ■ • , then  we  can  think  about  U operating  on  E(Xi) 

as  being  represented  by  an  infinite  matrix  G^.  This  matrix  is  referred  to  as  the  monodromy 
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matrix.  Its  eigenvalues  are  called  the  Floquet  or  characteristic  multipliers.  The  periodic 
solution  x(t)  of  (1)  is  stable  if  all  of  the  eigenvalues  of  U are  within  the  unit  circle  and 
unstable  if  there  is  at  least  one  with  a positive  real  part.  We  note  that  if  x(t)  is  an  exact 
periodic  solution  of  (1)  then  one  of  the  characteristic  multipliers  is  exactly  one. 

9.  Estimating  Characteristic  Multipliers 

In  this  section  we  assume  that  the  variational  equation  with  respect  to  the  approximate 
solution,  x(t),  can  be  written  in  the  form 

(66)  z(t)  — A{t)z(t)  + B(t)z{t  — u) 

where  A(t)  = A(t  + 2ir ),B(t)  = B(t  + 2n)  and  we  have  reintroduced  t to  make  the  oper- 
ator definitions  more  transparent.  Let  Z(t,s ) be  the  solution  of  (66)  such  that  Z(s,s)  — 
In,Z(t,s ) = 0 for  t < s where  In  is  the  n2  identity  matrix  on  Cn.  The  solution  Z(t,s) 
is  sometimes  referred  to  as  the  “Fundamental  Solution”.  Using  the  variation  of  constants 
formula  for  (66),  Halanay  [14]  shows  that  the  solution  of  (66)  for  the  initial  function  <j)  £ Ccj 
is  given  by 

(67)  z(t)  = Z(t.  0)^(0)  + j Z(t,a  + uj)B(a  + u))4>{a)  da 

J — Cj 

Define  the  operator 

(68)  (U(f))(s)  = z(s  + 2tt) 

where  4>  € C s £ [— u>,  0].  If  there  is  a non-trivial  solution  z{t)  of  (66)  such  that  z(t  + 2tt)  = 
pz(t)  then  p is  a characteristic  multiplier  of  (66).  If  we  combine  (67)  with  (68)  and  note 
that  z(a ) = <fi  for  a £ [ — a>,  0] , then  characteristic  multipliers  are  the  eigenvalues  of 

(69)  (Uct>)(s)  = Z(s  + 27t,  0)0(0)  + f Z(s  + 2tt,  a + 0)B(a  + Cj)<t>{a)  da 

J — tJ 

where  <f>  £ Cq.  Halanay  [14]  shows  that  we  can  restrict  s £ [ — u>,  0] . This  operator  is 
sometimes  referred  to  as  the  Monodromy  Operator. 

9.1.  Approximating  the  Fundamental  Solution  by  Spectral  Collocation.  In  this 
section  we  will  use  spectral  methods  to  compute  the  fundamental  solution  of  the  hnear 
homogeneous  delay  differential  equation  (66).  These  methods  are  well  known  for  collocating 
solutions  to  partial  differential  equations  and  boundary  value  problems.  See,  for  example, 
Gottheb  [11]  and  Gottlieb  and  Turkel  [12].  They  are  not  as  well  known  in  delay  differential 
equations.  In  this  section  we  use  a spectral  method  suggested  by  Bueler  [3]  and  Trefethen 
[31].  The  method  has  been  reported  earher  in  Gilsinn  and  Potra  [10]. 

The  computation  of  the  fundamental  matrix  used  in  the  monodromy  operator  (69)  re- 
quires the  computation  of  a solution  z(t)  of  (66)  on  some  interval  [a,  b\.  This  will  be  done  in 
a stepwise  manner.  We  first  find  a positive  integer  q such  that  a + qu)  > b.  Then  we  solve, 
at  the  first  step,  t £ [a,  a + u5], 

(70)  z\  (t)  = A(t)zi(t)  + B(t)zj(t  — u ;), 

where  z\(t  — u)  — 6(s)  for  some  function  (j)  £ Cb(a)  and  s = t — Q.  Thus  the  initial  problem 
becomes  an  ordinary  differential  equation.  Then,  on  [a  + uj,  a + 20}  we  solve 

(71)  z2(t)  = A{t)z2{t)  + B(t)z2{t  - 0), 

where  z2 {a  + 2)  = zi(a+£D),  z2{t  — u))  = zi(s)  for  s £ [a,  a + w],  s = t — Q.  Again  we  solve  (5) 
as  an  ordinary  differential  equation.  The  process  is  continued  so  that  on  [a+  (i  — l)a5,  a + iQ], 
for  i = 1, 2,  ■ ■ ■ . q, 

(72)  Zi(t)  = A(t)zi(t ) + B(t)zi(t  - 0), 

with  Zj(a  + (i  — 1)U)  = z,;_i(a+  (i  — l)cD).  We  then  define  z(t)  on  [a,  6]  as  the  concatenation 
of  Zi(t)  for  t £ [a  + (i  — 1 )Q,  a + iO]  and  i — 1,  2,  • • • , q. 
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Since  we  wish  to  use  a Cliebyshev  collocation  method,  we  will  shift  each  interval  [a  + (i  — 
1)2,  a + iu )]  to  the  interval  [— 1, 1],  For  t G [a  + (i  — 1)2,  a + iu],  for  i — 1,  2,  • • • , q,  we  have 
z G [—1,1]  provided 


(73) 


2_  (2a  + (2i  - 1)2) 

2 2 


For  z G [—1, 1]  we  have  t G [a  + {i  — 1)2,  a + iu]  provided 

2 (2a  + (2* -1)2) 

( < 4) 


U 

t = 1Z 


We  note  that  the  point  t G [a  + (i  — 1)2,  a + iu]  and  t — 2 G [a  + (i  — 2)2,  a + (i  — 1)2]  are 
translated  to  the  same  z G [—1,1].  This  is  clear  from 


(75) 


2 ^ (2a  -f  (2z  — 3 )u  2 (2a  + ( 2i  — 1 )u 

~{t  -u) = —t 

u u to  u 


Therefore  we  can  shift  the  iterated  delay  problems 

(76)  Zi(t ) = A(t)zi(t)  + B(t)zi(t  - 2), 

for  t G [a  + (i  — 1)2.  a + iu]  and  i = 1,2,---  , q,  into  iterated  ordinary  differential  equations 

(77) 

where,  for  t G [a  + (i  — 1)2,  a + iu]  and  associated  z G [—1.  1], 

Ui(~  1)  = Ui- i(l) 

Ui(z)  = Zi(t) 


u'iiz)  = -Ai(z)ui(z)  + -Bi(z)ui-i(z) 


(78) 


Mz) 

Bi(z) 

Ui-l(z) 


A(t ) 

m 

Zi(t  — 2) 


■uo(^)  = zi(t  - 2)  = <j>{t  - 2) 


The  initial  function  is 

(79) 

for  t — 2 G [a  — 2,  a]. 

We  can  now  approximate  the  fundamental  solution  for  (66)  on  [a.  6]  by  first  solving  the 
iterated  differential  equations  (76)  subject  to 

Ui{- 1)  = Ui-l(l) 

(80)  u0(z)  = 0,  z E [—1,  l] 

«l(-l)  = In 

where  In  is  the  n x n identity  matrix.  We  follow  the  spectral  method  given  in  Bueler 
[3]  in  that  the  fundamental  solution  is  solved  for  in  n passes  of  the  iteration  process  with 
«i(  — 1)  = ej.  where  ej  — (0,  ■ ■ • , 1,  • • ■ , O)7  with  1 in  the  j th  element,  j = 1,  2,  • • • , n. 

To  begin  the  solution  process  we  take,  for  some  positive  integer  N.  the  Chebyshev  points 


(81) 


kn 

Vk  = cos  [ — 


on  [—1, 1],  for  k = 0, 1,  ■ • ■ , N.  The  benefits  of  using  these  points  has  been  discussed  by 
Salzer  [24].  The  Lagrange  interpolation  polynomials  at  these  points  are  given  by 

N 


(82) 


tw = n 


Vk 


k^j 


We  have  lj{rjk)  — djk-  Then  on  [—1, 1]  we  set 


N 


ui(z)  = YlUi  fa)  h 

3=0 


(83) 
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We  also  need  to  form 


TV 

(84)  u'^z)  = YUi  (rjj)  Z'-(z) 

j= o 

At  the  Chebyshev  points  we  will  designate 

(85)  Dkj  = lj(rjk) 

The  values  for  these  derivatives  are  given  in  Gottlieb  and  Turkel  [12]  or  Trefethen  [31]  but 
we  state  the  values  for  D here  for  completeness.  The  derivations  are  given  in  Section  15, 
Appendix  1. 

2N2  + 1 


1! 

O 

o 

Cl 

6 

Dnn  = 

—Doo 

(86) 

Dn  = 

-Vj 

1,2,---  , N 

2(1  - 

2\  ’ J 

Vj) 

Dij 

Ci{- 1 

y+j 

Cjim  ~ 

-Vj) 

o 

II 

•<s> 

"tk 

■ , N where 

(87) 

Ci 

-ft 

i = 0 or  TV; 

otherwise. 

For  notation,  let 

Ui{z) 

= ( 

\T 

j V'in)  i 

(88) 

Mz) 

= 

p.q=  I.-*-  .n 

B,(z) 

= 

l p.g=l,---  ,n 

We  then  write  the  collocation  polynomial  of  Ujr,r 

= 1,-  - ■ ,n, 

TV 


(89)  uir{z)  = Y,  u-ttk(z) 

at  the  Chebyshev  points  (81)  to  get 


k= o 


(90) 


uir  (‘Vj  ) 

= w 

rj 

N 

Uir(Vj ) 

= x> 

fc=0 

Hj  — l.r 

= w(ir 

rj 

,(*) 


yjk 


The  initial  conditions  for  the  iterated  differential  equations  are 


(91) 


Hir  (Vn)  — V-i  — l.r  iVo)  i 


or 

(92) 


,0) 

■rTV 


= W 


0-1) 

r0 


for  r = 1,  • • • , n. 


The  discretized  differential  equations  are  then  given  by 
(93) 
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(94) 

Wj  = 

( (i) 

\wio  • • • 

(i)  (i) 

W1NW20 

(i) 

' ‘ ' W2N  ■ 

( i ) 

' ' Wn0 

( (i— 1) 

(i- 1) 

0-1) 

(i- 1 

Wj- 1 = 

(^10 

• ' ' W1N 

W20 

' W2N 

■ W. 


;nN 


for  j = 0, 1,  ■ • • , N — 1.  These  provide  nN  equations  but  n(N  — 1)  unknowns.  The  other  n 
equations  come  from  the  initial  conditions.  We  define  the  following  vectors 

W \T 

nN  J 

r/ir1) 

“ 77.0 

Then  we  can  write  the  iterated  differential  equation  as 

(95) 

where  D = D ® In , the  Kronecker  product,  and  each  D is  given  by 

Z?oo  • • • Don 

(96)  D = 


~ id  ~ id  ~ 

Dwi  = —AiWi  + —BiWi 


Dn- i.o 
0 


Dn-i.n 

1 


The  unit  in  the  lower  right  introduces  the  initial  condition,  wr^,r  = 1,  - , n,  equation. 

Thus  D is  formed  by  n blocks  of  D down  the  diagonal. 

The  matrix  Aj  is  given  by 
(97) 


Ai  = 


1 

o 

0 

0 • 

■ 0 • 

• 0 

0 

0 • 

• o' 

0 

0 

0 ■ 

• 0 ■ 

• 0 

0 

0 ■ 

• 0 

0 

0 

^ii(ihv-i) 

0 • 

■ 0 • 

• 0 

0 

A-iniVN-l) 

0 ■ 

• 0 

0 

0 

0 • 

• 0 • 

• 0 

0 

0 • 

• 0 

AnliVo) 

0 

0 • 

■ 0 • 

• 0 

A(nn{r)o) 

0 

0 • 

• 0 

0 

0 

0 • 

• 0 • 

■ 0 

0 

0 • 

• 0 

0 

0 

^nl(^-l) 

0 • 

• 0 • 

■ 0 

0 

Ann(^)N—  1 ) 

0 • 

• 0 

0 

0 

0 • 

• 0 ■ 

• 0 

0 

0 • 

• 0_ 

B,  is  structured  in  a similar  manner  except  every  (N  + l)th  row  includes  an  element  2/Q 

to  take  care  of  the  initial  condition.  Thus 

(98) 


B,  = 


'b[  %o) 

O' 

0 ■ 

• 0 ■ 

• 0 

0 

0 • 

• o' 

0 

0 

0 • 

• 0 ■ 

• 0 

0 

0 • 

• 0 

0 

0 

bS(vn- i) 

0 • 

■ 0 • 

• 0 

0 

B$(vn- i) 

0 • 

• 0 

2_ 

p 

0 

0 • 

■ 0 • 

• 0 

0 

0 • 

• 0 

bS(vo) 

0 

0 • 

• 0 ■ 

• 0 

B{n% 70 ) 

0 

0 ■ 

• 0 

0 

0 

0 • 

• 0 ■ 

• 0 

0 

0 • 

• 0 

0 

0 

B%(m-i) 

0 • 

■ 0 • 

• 0 

0 

Bnn(jjN  — 1 ) 

0 ■ 

• 0 

0 

0 

0 • 

• 0 • 

• 0 

2 

u) 

0 

0 ■ 

• 0_ 

The  linear  equation  (95)  can  be  solved  for  Wj  by  setting 

-l  ~ 


(99) 

and 

(100) 

for  i = 2,3,- 


Mi  = [D  - — Aj  -B, 


Wj  = MjWi— i 
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To  solve  for  w\  for  the  fundamental  solution  we  need  to  solve 

(101)  u[(z)  = ^A1(z)u1(z) 
for  2 G ( — 1. 1]  and 

(102)  Ul(-l)=/n 

That  is,  we  solve  n problems  at  each  iteration,  one  for  each  of  the  initial  conditions  el5  where 
e,  is  the  standard  basis  vector  with  a unit  in  the  xth  element  and  zero  elsewhere.  For  the 
moment  we  set  the  initial  vector  as 

(103)  w0  = (0  • ■ • noiO  • ■ • 1x020  ■ • ■ u0n)T 

where  uor,r  = 1,  • ■ • ,n,  is  placed  in  each  of  the  ( N + l)th  elements  and  zero  elsewhere. 
Then  from  the  previous  construction  of  D and  Ai  we  have 

(104) 

Given  that  we  have  computed 

N 

(105)  Uir{z ) = ^ Wrklk(z ) 

fc= 0 

on  [—1,1]  for  r = 1,  • • • , n we  can  compute  the  result  for  t G [a  + (x  — 1 )h,  a + ih\  by  setting 

(106)  zir(t)  = uir(z) 


for  r = 1,  • • • , n,  where 


(107) 
or 

(108) 


2 (2a  + (2x  - 1)£) 

id  id 


•(*)  =J2Wrklk(~t  - 

— id 

k= 0 


2 (2a  + (2i-l)D) 


The  initial  condition  is 


(109)  uir(r]N)  = Ui-ltr(r)o). 

But  on  [a  + (x  — l)cD,  a + iQ],  zjv  = —1  corresponding  to  t — a + (x  — l)u)  and  on  [a  + (x  — 
2 )u;,a  + (x  — 1)lD],  zq  — 1 corresponding  to  t = a + (x  - l)u),  so  that 

(110)  Zir  (aT  (x  — 1 )ca)  — — i,r  (^  T (x  — 1 )tn) 

9.2.  Estimating  Monodromy  Operator  Eigenvalues.  To  approximate  the  monodromy 
operator  (69)  we  will  require  a quadrature  rule  that  satisfies 


(111) 


P+ 1 pO 

^2  Vkf  (sk)  / /(s) 


as  P ->  oo  for  each  continuous  function  / £ Cj,.  The  rule  is  satisfied  if 


P+i 

(112)  M < Q, 

fc= l 

for  some  Q > 0 and  P = 1,  2,  ■ • • . This  is  satisfied  by,  for  example,  Trapezoidal  or  Simpson 
rules. 

Let  — Cj  = s\  < S2  < ■ ■ ■ < = 0,  and  define 


p+i 

(113)  (U (fr)  (s)  = Z (s  + 2k . 0)0(0)+  y v^Z  (s  + 27r,Sfc  + u>)  B (sk  -\-  w)  <t>  (sd 

k=i 


for  4>  G Cqj. 
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Then,  for  each  Sj  € [— u>,0], 


p+i 


(114)  (170)  ( Si ) = Z(sj  + 271, 0)0(0)  + ^ WjZ(sj  + , s j + Cj)B(s j + cj)<fi(sj) 

i=i 

Since  sp+i  = 0,  (114)  can  be  rewritten  as 

p 

(U</>)  (si)  = E WjZ(s,  + 271,  Sj  + u})B(sj  + u)<j>(sj ) 

j=i 

(115)  + (Z(sj.  + 2n,0)  + wp+\Z(si  + 2tt,uj)B(uj))  <f)(sp+i) , 

where  Z(s,a)  is  the  fundamental  matrix  of  (66).  Equation  (115)  can  be  put  in  matrix  form 


/ m)(S j \ 

- U\.\ 

Uij  ••• 

U\.p+\ 

(116) 

(U4>)(si) 

= 

Ui,  1 • • • 

Uij 

Ui,P+1 

{(U4>)(sP+  l)J 

.C^p+i,i  ■ • ■ 

Up+i.j  ■ ■ ■ 

Up+i,P+i_ 

where  the  block  elements  for  i = 1,  ■ • ■ , P + 1,  j = 1,  ■ • • , P are  Uij  = wjZ(s,  + 271,  s7  + 
uj)B{sj  + Cj).  The  block  elements  in  the  last  column  of  the  matrix  are  given  by  C/j.p+i  = 
Z(sj  + 27r,  Q)  + wp+\ZiySi  + 2'K ,uj)B(Cj)  for  i — 1,  • • ■ , P+  1.  The  relevant  eigenvalue  problem 
becomes 


‘ Uip  ■■■ 

Upj  ••• 

Ui,p+1 

( (0)(s l)  ^ 

(117) 

Ui.  1 • • • 

Uij 

Ui.p+i 

= A 

WiSi) 

Up+i.i  • • • 

Up+i,j  ■■■ 

Up+\.p+\_ 

V(0)(sp+ 1)/ 

10.  Determining  Solutions  of  the  Adjoint  Equation  Associated  with  Multipliers 

of  the  Variational  Equation 

In  order  to  estimate  a in  (33),  let  t & [0,  27i]  and  ip  be  the  initial  function  defined  on 
[271,271  T ui] . The  formal  adjoint  equation  from  (13)  is  given  by 

(118)  v0 (t)  = -v0(t)A(t)  - v0(t  + Q)B{t  + Q), 

where  vq (t)  is  a row  vector.  Ordinarily  solving  the  adjoint  equation  would  require  a ’’back- 
ward” integration.  However,  Halanay  [14]  showed  that  the  solution  of  the  formal  adjoint 

(118)  on  [0,  27i]  is  given  in  row  vector  form  by 

(119)  vo(t)  = ip(2n)Z(2n,t)  + f ip(a)B(a)Z(a.  — Q,t)  da. 

■h-K 

The  significance  of  this  representation  is  that  only  a “forward”  integration  is  required  to 
solve  for  the  fundamental  solution  Z.  Let  0(s)  be  a continuous  row  vector  function  defined 
on  [— u5,  0].  Halanay  [14]  then  defined  the  operator 

(120)  (Ucp)  (s)  = (p(—u))Z(2n,s +Q)  + I (p(a)B(a  + Q)Z(2tt  + a,  s + Q)  da. 

J —Q 

Note  the  relationship  to  the  monodromy  operator  (??).  He  also  gave  an  associated  operator 
V,  defined  on  [2n,  2n  + u5],  as 

(121)  (V'lp'j  (s)  = y(s  — 2n,'ip)  = ip(2Tr)Z(2TT,  s — 2tt)  + J tp(a)B(a)Z(a— 0,  s — 2n)  da. 

Again,  note  the  relationship  to  (119).  He  further  showed  that  an  eigenvalue  po  of  V is 
associated  with  a 1 / p0  multiplier  of  the  formal  adjoint  equation  (118),  the  eigenvalues  of 
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U,  U,  V are  all  the  same,  and  the  eigenvectors  of  U,  V are  related  by  0(s)  = if)(s+2ir+Q),  s G 
[w,0].  Although  V is  the  operator  associated  with  (119),  the  fact  that  the  eigenvalues  and 
eigenvectors  of  V and  U are  the  same  allows  algorithms  developed  for  the  characteristic 
multipliers  in  Section  9 to  be  easily  modified  to  compute  the  eigenvalues  and  eigenvectors 
for  (120).  In  particular  we  again  partition  the  interval  [— tD,  0]  into  P equal  intervals  of 
length  A = u)/F  by 

(122)  —1 d = s\  < S2  < ■ ■ ■ < sp+i  = 0 

We  use  a Simpson  integration  method  to  write,  with  the  same  weights  as  previously, 

(s)  = 0 (si)  (z(27r,  s + lD)  + w\B  (si  + Q)  Z (2tt  + si,  s + 2)^ 

p+i 

(123)  + 'y  ^ (ft  (sj)  (wtjB  {sj  + lo)  Z (27t  + Sj,  s + ui)^ 

1=2 

To  solve  the  adjoint  equation  in  row  form  on  [0,  27t],  we  need  only  compute  the  eigenvector 
of  U associated  with  the  characteristic  multiplier  of  U . The  eigenvector  is  then  substituted 
into  equation  (119).  From  the  previous  section  there  are  likely  to  be  two  complex  conjugate 
eigenvalues,  po  and  p0,  associated  with  complex  conjugate  eigenvectors  vq  and  ~o  of  (119), 
so  by  linearity  of  (118)  the  real  part  forms  a discretized  solution  of  (118)  and  thus  a single 
real  independent  solution.  We  will  continue  to  call  the  real  part  of  this  eigenvector  0 so  that 
for  t € [0,  27t] 


v0{t)  = cj)(si)  (z(2tt,  t)  T W\B  (si  T cj)  Z (si  T 27t,  t)^ 

P+i 

(124)  + y ^ <p  ( Sj ) [wjB  ( Sj  + uj)  Z ( Sj  + 27t,  t)^j 

3= 2 


The  j-th  block  column  of  170  is  given  by 

(125) 


( Sj ) = 0(Sl),-"  ,0(S;),--  - ,0(Sp+l) 


Z(2tt,  Sj  T id)  T w\B)s\  u)A  (3^  + 2i r,  s j T u^) 


W i B ( .S { 4~  Ld)Z(Si  T 27T,  Sj  T id) 


wpjr\B{s]\[j-i  + cd)Z (s 7v+i  T 27 r,  Sj  + id) 


The  eigenvector  0 of  the  matrix  on  the  right  associated  with  the  multiplier  of  the  variational 
equation  is  computed  and  substituted  into  equation  (124)  to  give  the  value  of  vq (t)  on  [0.  2tt\. 

To  compute  a we  need  to  estimate  (13).  Again  we  use  a Simpson  rule.  We  partition 
[0,  27 r]  by  equidistant  intervals  0 = t\  < < ■ • ■ < fo+i  = 27r  and  set  h = 2i t/O,  where  O 
is  an  even  integer.  Again  the  weights  will  be  set  as  u\  = uo+i  = h/ 3,  otherwise  Uk  = 4h/3 
if  k is  even,  Uk  = 2/i/3  if  k is  odd.  We  then  set 


(126) 


a = 


(o+i 

y,  ukv0  (tk)  j 

\k=i 


-i 


where  Vq  is  a row  vector  and  J is  a column  vector. 
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Using  the  same  partition  of  [0.  27t]  we  can  compute  vq (t)  at  each  tk  as 

(127) 

Z(2n,tk)  + w\B(s\  + 27r  + Q)Z(si  + 2tt,  tk) 

Vo (tk)  = |0(«i),  • • • ■ ■ ■ ,4>(sp+i)}  uijB(si  + 27T  + Q )Z(si  + 2n,tk) 

wp+\B{spjri  + 27t  4-  Q)Z(sp+i  + 27t,  tk) 

We  finally  normalize  vq  so  that  |uo|2  = 1- 

11.  Estimating  the  M Parameter 

From  Halanay  [14]  the  variation  of  constants  formula  for 

(128)  z(t)  = A(t)z(t ) + B(t)z(t  - Cj)  + /(f), 
where  t £ [0,27r],  is  given  by 

(129)  z(t)  = Z(t,  0)0(0)  + f Z(t,a  + u>)B(a  + t2>)z(a)  da  + I Z(t,a)f(a)da. 

J -Cj  Jo 

The  27t  periodic  initial  function  condition  with  s € [—<£>,  0]  is 
0(s)  = Z(s  + 27t,  0)0(0) 

/•0  rS-\-  27T 

(130)  + / Z(s  + 27t,  a + u>)B(a  + Cj)<j)(a)  da  + / Z(s  + 2n,a)f(a)  da. 

J —Cj  J 0 

The  first  step  in  computing  M involves  relating  0 to  /.  Let  |0|  = sup_ti<s<0  |0(s)|  and 
similarly  for  |/|  on  [0,2tt].  To  eliminate  0(0)  from  (130),  set  s = 0 in  (130)  and  solve  for 
0(0)  as 

r° 

0(0)  = / (/  — Z(27t,0))  1 Z(27t,  a + Ci)B{a  + u;)0(a)  da 

J — Cj 


o2tt 


(131) 


+ / {I  ~ Z{2n,0))-iZ(27r,a)f{a)da 

Jo 

Substitute  (131)  into  (130)  and  combine  terms  as 
r0 


(132) 


0(s)  = f [Z(s  + 2n,  0)(7  — Z(2n.  0))~1  Z(2n,  a + Cj) 

J —Cj 

+ Z(s  + 27t,  a + w)]  B(a  + u))0(ct)  da 

p2n 

+ / [Z(s  + 2tt.0)(J  - Z(27T,0))-1Z(27T,a)  + Z(s  + 27r,a)]  f(a)da. 
Jo 


where  s £ [ — u>,  0] . 

Let  —Cj  = s\  < S2  < ■ ■ ■ < sp+ 1 = 0,  ds  = *p,  and  0 = fi  < t2  < ■ • • < to+i  = 2tt,  df  = 


jf.  We  can  discretize  (132)  by  setting 


(133) 

( 

where 

(134) 

P+l  0+1 

0(*i)  = X!  ) 

i=i  fc=i 


TZ  (3;  T 27t,  Sj  T ca)  j L?  (3j  -)-  uj) 
ua-  (sj  + 27t,  0)  (I  — Z(27t.  0))_1 
+Z  ( Si  + 2tt,  tk)  j 
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In  vector  matrix  form  (133)  can  be  written 


( 0(Si)  \ 

( 4>{si)  \ 

/ fit  1)  \ 

. 

+ hA  : 

\4>{8P+1)j 

\<t>{sp+ 1)/ 

\f(to+i)J 

(135) 

Using  a generalized  inverse  we  can  solve  for  the  (j>  vector  with  minimum  norm  by 

(136) 


f 4>{sx)  ^ 

= {I-H1)+H2 

( fit  1)  \ 

\<t>{sp+i)/ 

U(io+i)/ 

In  the  second  step  the  value  of  </>(0),  given  by  equation  (131),  is  substituted  into  equation 
(129)  and  terms  combined  to  give 


:{t) 


137) 


/ [Z(t,0)(I-  Z(2^,0))-1Z(27r,a  + th) 

J — Cj 

+ Z(t,a  + tu)]  B(a  + u j)</>(a)  da 


+ 


/»27T 

/ [Z(t,0)(I-  Z(27T,0))-1Z(27r,a)  + Z(t,a)]  f{a)da. 

Jo 


This  can  be  discretized  by  setting 

P+ 1 


O+i 


(138) 


(tfc)=  J]  H3(k,i)4>(Si)+  ]T  H4(k,j)f(tk) 

J=1 


2—1 


where 


H3(k,i) 

= v,  [z  (tk,  0)  (7  - Z(2tt,  0))_1Z  (2? r,  s,  + w 

139) 

+ Z (tfc,  S,  + th)  j 77  (s,'  + ii>) 

= Uj  [z  (tfc,  0)  (7  - Z(2tt,  0))-jZ  (2tt,  fj) 

+Z  ( tk , tj)  j 

In  vector  matrix  form  (138)  can  be  written 


/ ^(*1)  \ 

( <p{s\ ) \ 

( fit  1)  \ 

(140)  : 

= 77.3 

+ 774 

\z(to+ 1)) 

V0(«p+i)y 

\f{to+i)J 

By  substituting  (136)  into  (140)  we  have 


(141) 


/ z{h)  \ 
\z(to+i)j 


[h3  (I  - Hi)+  H2  + 774] 

/ /(7i)  \ 

U(7o+i)/ 

Therefore 


(142)  \z\  < M\f\, 

where  M — \\h3  (I  — H\)+  H2  + 774 II  . 
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Time 


Figure  1.  Residual  Error  of  Approximate  Solution  for  the  Van  der  Pol  Equation. 


12.  Application  to  a Van  der  Pol  Equation  with  Delay 

In  this  section  we  will  apply  the  main  theorem  to  approximate  the  limit  cycle  of  the  Van 
der  Pol  equation  with  unit  delay,  given  by 

(143)  x + A (x(t  — l)2  — l)  x(t  — 1)  + x = 0. 

Since  the  period  of  the  limit  cycle  is  unknown  we  introduce  an  unknown  frequency  by 
substituting  t/uo  for  t to  obtain 

(144)  ui2x  + u>\  (x(t  — co)2  — l)  x(t  — u>)  + x = 0, 

for  t E [0.  27t] . To  compare  with  an  approximation  result  obtained  for  ordinary  differential 
equations  in  Stokes  [27].  we  take  A = 0.1. 

The  first  step  was  to  estimate  an  approximate  2yr-periodic  solution,  frequency  and  residual 
to  (144).  By  using  Galerkin’s  method  described  in  Section  7 the  following  approximate 
solution  was  obtained 

x(t)  = 2.0185  cos(i) 

+ 2.5771  x 10-3  sin(2t)  + 2.5655  x 10“2  cos(2 1) 

(145)  + 1.0667  x 10-4  sin(3i)  - 5.2531  x 10“4  cos(3 1) 

- 7.1780  x 10“6  sin(4t)  - 2.2043  x 10'6  cos(4t), 

Co  = 1.0012. 

where  we  have  displayed  only  the  first  few  harmonics.  This  solution  was  estimated  based  on 
11  harmonics,  40.000  sampled  points  over  [0,  27r],  and  100  Chebyshev  extreme  points  (81). 
The  residual  was  estimated  by  substituting  (Co,x)  from  equation  (145)  into  equation  (144) 
and  finding  the  maximum  of  the  absolute  values  of  the  residuals  obtained  in  the  interval 
[0.27t],  The  result  was  r = 3.1086  x 10~15.  This  residual  is  significantly  better  than  the 
one  given  in  Stokes  [27].  The  distribution  of  the  residuals  for  the  current  case  is  shown  in 
Figure  1.  The  phase  plot  of  the  approximate  solution  is  shown  in  Figure  2.  For  t E [0.  2 tt] 
we  can  then  hnmediately  estimate  |x|  < 2.0436,  |x|  < 2.0279.  |x|  < 2.1165. 

In  the  second  step,  the  values  of  the  constants  B and  1C  were  obtained  in  a straightforward 
manner  from  the  variational  equation  about  the  approximate  frequency  and  solution  given 
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Phase  Plot  of  Approximate  Solution 


X 


Figure  2.  Phase  Plot  of  Approximate  Solution  for  the  Van  der  Pol  Equation. 


by 

(146)  z{t)  = A(t)z(t)  + B(t)Z(t  — lo), 


where 


zi 


A(t)  = 


0 

— 1/ch2 


B(t)  = 


0 0 
— 2(A/J;)xi(t  — Cj)x2{t  — u)  (A/T)  (l  — x\ (t  — u)2) 

We  use  the  fact  that  the  natural  norm  of  a matrix,  H , associated  with  a vector  norm 


\x\  = maxi<i<„  \Xi\ 

show  that 


is 


\H\  = maxi<j<n^)”=1  \hij\.  With  this  definition  it  is  not  hard  to 


(147) 


\dX(x-,4>)\  < 

0 1 

— l/u>2  — 2(A/ut)5'l(f  — Cu)X2(t  — u)  (A/th)  (l  — fi(t  — th)2) 
< 2.3776101. 


\4>\, 


Therefore,  for  A = 0.1.  B = 2.3776.  Working  conservatively  within  the  domain  D — 
{x  £ C[0,  27t]  : \x  — x\  < 1}  it  is  not  hard  to  show  that 

\dX(xcj  + V’ i;  <Ai)  - dX{x^  + lp2',4>w)\ 

(148)  < (6A/c2>)  (1  + |£|)  |0i  — 02|  |0|- 

Then  from  (145)  and  (148)  we  can  estimate  1C  — 1.8157  and,  from  (23),  we  can  estimate 
J(f,th)j  < 2.7546. 

Next,  we  can  estimate  the  characteristic  multipliers  of  the  variational  equation  relative 
to  the  function  x(t).  For  the  quadrature  steps  in  Sections  9 and  10  P and  O were  taken 
as  200  and  1200  respectively.  These  gave  mesh  widths  of  about  1/200  on  both  [— u),0]  and 
[0.  27r].  Using  the  method  of  Section  9 we  computed  two  simple  conjugate  eigenvalues  with 
magnitude  1.0430.  All  of  the  other  eigenvalues  have  magnitudes  near  zero.  These  are,  of 
course,  the  eigenvalues  of  the  monodromy  operator  U.  The  fundamental  matrix  Z in  (69) 
is  computed  using  the  collocation  method  of  Section  9.1  (See  Figure  3).  The  monodromy 
operator  is  formulated  as  in  Section  9.  The  eigenvalues  of  the  monodromy  operator  U are 
plotted  in  Figure  4.  Note  that  the  significant  complex  conjugate  eigenvalues  are  near  the 
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Z(1,1) 


Z(2.1) 


tt 

Z(1,2) 


tt 

Z(2,2) 


tt  tt 


Figure  3.  Fundamental  Matrix  for  the  Variational  Equational  relative  to 
the  Approximate  Solution  for  the  Van  der  Pol  Equation. 


Eigenvalues  for  U 


Figure  4.  Eigenvalues  for  the  Monodromy  Operator 


unit  circle  but  are  not  exactly  on  it.  This  is  due  to  the  fact  that  (145)  is  only  an  approximate 
solution.  The  eigenvalues  are  complex  conjugates  because  the  left  hand  matrix  in  (117)  is 
real  and  non-symmetric  since  the  fundamental  solution  Z is  non-symmetric  (See  Figure 
3).  We  can  confirm  that  the  eigenvalues  of  the  operator  U are  the  same  as  those  of  U . 
Graphically  this  is  shown  in  Figure  5. 

In  the  next  step  we  estimate  the  parameter  a using  the  methods  of  Section  10.  The 
solution  of  the  adjoint  to  the  variational  equation  was  computed  using  equation  (127)  and  the 
parameter  a in  (33)  was  estimated  by  simple  quadrature,  with  A = 2-n/O  for  a sufficiently 
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Eigenvalues  of  U— tilde 


Figure  5.  Eigenvalues  for  U 


large  mesh,  0 = t\  < t2  < ■ ■ ■ < to- 1 = 27t,  as 


(149) 


a = 


A 


o+i 

y{ti)J{x,u)(ti) 


-i 


f=i 

The  absolute  value  of  a is  estimated  as  3.3547. 

If  we  now  apply  the  methods  of  Section  11,  using  A(t)  and  B(t)  defined  in  equation  (146), 
we  can  estimate  M = 2.7618  x 102.  These  results  allow  us  to  estimate  Ao,  Ai  and  A2  in 
Lemma  4.3  as  Ao  = 8.4091.  Aj  = 6.6736  x 103,  and  A2  = 3.1720  x 104.  Note  the  magnitude 
of  the  parameters. 

With  the  estimates  above  we  can  compute  Fi  — 2.5941  x 109.  F2  = 1.0798  x 1010  from 
(202)  and  (205)  respectively.  Then  we  compute  8 = 4.6305  x 10-11  from  (54).  Then 
F1S2  = 5.5623  x 10~12  is  less  than  8 and  F28  — 0.5.  Furthermore  r < 8.  Therefore,  the 
conditions  of  the  main  theorem  are  satisfied  and  we  can  conclude  from  Theorem  6.1  that 
there  exists  an  exact  solution  x*  and  an  exact  frequency  uj*  of  equation  (144)  such  that 
\x*  - x\  < 1.2361  x 10"6  and  |cu*  - u\  < 7.7877  x 10“10. 


13.  Conclusions 

Although  there  seem  to  be  a large  number  of  parameters  to  be  computed  and  inequalities 
to  be  tested  in  order  to  produce  the  final  error  estimates  the  process  is  feasible.  All  of  the 
steps  can  be  completed  within  a single  code.  The  current  code  in  Appendix  3,  Section  18,  has 
also  been  built  around  the  example  in  Section  12  and  would  have  to  be  generalized  for  other 
applications,  but  the  code  provides  a template  on  which  to  proceed.  From  the  computational 
point  of  view  the  longest  compute  times  involve  the  construction  of  the  block  matrices 
(116)  and  (125).  Computing  the  approximate  solution  and  the  fundamental  solution  of  the 
variational  equation  is  relatively  fast  compared  to  these  matrix  constructions.  It  behooves 
anyone  wishing  to  apply  the  methods  of  this  paper  to  spend  some  effort  vectorizing  the 
matrix  construction  algorithms  in  Sections  9.1  and  10  as  much  as  possible. 

The  parameter  M in  the  Fredholm  Lemma  3.2  is  a significant  parameter.  From  the 
example  above,  it  is  clear  that  it  would  be  desirable  to  obtain  as  small  a value  for  M as 
possible,  since  its  magnitude  affects  the  A;,  i = 1,2  parameters  and  Ai  appears  in  the 
final  error  estimates.  I11  particular,  in  the  example  above,  the  effect  of  M causes  a very 
fine  residual  r for  the  approximate  solution  (145)  to  produce  a pessimistic  error  estimate 
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between  the  approximate  solution  and  the  exact  solution  in  the  end.  From  (40)  the  critical 
parameter  Ai  is  linearly  dependent  on  M. 


14.  Disclaimer 

Certain  commercial  software  products  are  identified  in  this  paper  in  order  to  adequately 
specify  the  computational  procedures.  Such  identification  does  not  imply  recommendation 
or  endorsement  by  the  National  Institute  of  Standards  and  Technology  nor  does  it  imply 
that  the  software  products  identified  are  necessarily  the  best  available  for  the  purpose. 


15.  Appendix  1 

In  this  appendix  we  present  the  derivation  of  the  differentiation  matrix  (86).  The  deriva- 
tion is  based  on  a discussion  of  pseudospectral  Chebyshev  methods  given  in  Gottlieb  et  al. 
[13].  although  a full  derivation  of  the  differentiation  matrix  is  not  given. 


Lemma  15.1.  For  some  positive  integer  N let  the  Chebyshev  points  be  given  by 


(150) 


on  [—1, 1],  for  k — 0,1,---  , N.  The  Lagrange  interpolation  polynomials  at  these  points  are 
given  by 


(151) 


We  have  lj(r}k)  = $jk-  At  the  Chebyshev  points  designate 
(152)  Dkj  = I'jirjk) 

The  values  for  these  derivatives  are  then  given  as 

2 N2  + 1 


(153) 


D oo  = 
Dnn  = 
Djj  = 

Dij 

, N where 


6 

— Doo 
~'h 

2(1-^) 
c,(-l)i+J 
CjiVi  ~ Vj) 


,3  = 1,2, 


2,  i = 0 or  N; 
1.  otherwise. 


,N-  1 


fori  ^j,i,j  = 0. 

(154)  a = 

Proof:  The  Chebyshev  polynomial  of  degree  N is  given  by 


(155) 

for  z £ [-1,  1], 

Define  the  polynomial 

(156) 


Tn(z)  = cos  (fV cos  x) 


9j{z) 


(1-22)T^(2)(-1)1+1 

CjN2  ( 2 — Zj) 


for  j = 0,  • • • ,N  and  Co  = cn  = 2,  Cj  = 1 for  1 < j < N — 1.  Since  T'N  (zj)  will  be  shown 
below  to  equal  zero,  T'N(z)/  (z  — Zj)  is  a polynomial  of  degree  N — 2 so  g3  (z)  is  a polynomial 
of  degree  N.  Thus,  if  we  can  show  that  gj  (zk)  = 6jk  for  k = 0.  • • • , N,  then  by  uniqueness 
9j(z ) = 
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We  first  need  to  compute  the  following  derivatives. 

—N  sin  (TV  cos-1  z) 


T'n(z)  = 


vT^T2 


—N2(l  — z2)1/~cos(Ncos  1 z)  — Nz  sin  (IV  cos  1 z) 
(1  - z2)3/2 

T'^(z)  = —N2  sin  (iVcos-1  z)  Af  (l  — z2)~l/2  (l  - z2)-1 

+ cos  (TV  cos-1  z)  (-1)  (l  — z2)  “ (— 2z)J 
—N  [sin  (iVcos-1  2)  (l  — z2)  3/~ 

+ Z COS  [N  COS-1  z)  iV  (l  — 22)  1'"(l—  Z2)  3,/“ 


(157) 


9j(z)  = 


+z  sin  (iV  cos 


-1 


(1-z2)  5/2  ( — 2z) 


CjlV2 


(~2z)T'n(z)  | (1  ~ z2)  T&(z) 


Z Z q 


z — 


+ 


(i-*2)n 


V v 


(z-z,r 


(-i)j+1 

CjlV2 
AT  (-1  - 


./Vz  SHI  (Af  COS  1 2)  A"2  COS  (iV  COS  1 z) 


(z  - Z,-)  (1  - Z2)1/2 


sm  (A  cos 


(*"*) 


We  will  first  establish  that  c/j(z)  = lj(z).  Clearly,  since  cos-1  Zfc  = ktr/N , T/v  (zfc)  = 0, 
and  therefore,  for  fc  ^ j,  k ^ 0,  Ah  j 7^  0.  Ah  gj  (zfc)  = 0.  For  k - j.  j ^ 0.  N,  using  T'N(z) 
and  L ‘Hospital’s  rule  for 


(158) 


lim 


sin  (AT  cos  1 z)  A(  — 1)J 


(l-z2)1/2' 


we  have  gj  ( Zj ) = 1.  For  j = 0,  2o  = 1 so  that 

( — 1)J+2  (l  - z2)  l,'~  sin  [N  cos-1  z) 


(159) 


9o{z)  = 


2N(z  - 1) 

For  z = Zfc,  k 7^  0,  go  (. Xfc ) = 0.  Again  apply  L'Hospital’s  rule  to  show 

z sin  (AT  cos-1  2) 


(160)  go  (zo)  = lim 

27V  z->  1 


Ar  cos (A  cos 


(1-z2) 


2 11/2 


= 1. 


For  j = Ah  z/v  = — 1 and  (zfc)  = 0 for  k = 0,  1,  • ■ ■ , A — 1.  For  fc  = 0.  use  L Hospital's 
rule  to  show 


(161) 


gN  ( zn ) — 


-i)^2  .. 

— ;!“i 


: sin  (A  cos  1 z) 
(l-z2)1/2 


+ A cos  (A  cos  z) 


= 1. 


Therefore,  gj(z)  = Zj-(z). 

We  now  construct  the  entries  in  the  differentiation  matrix  (177).  These  are  given  by 
Djk  — g'k  (zj)  f°r  A k = 0,  1.  -■  , Ah  For  A:  7^  j,  k 7^  0,  A.  since  sin(fc7r)  = 0 and 
cos(kn)  = ( — l)fc, 


5'  (**)  = 


Cfc(-l)J  + 1 
cj  {zk  — Zj ) 


(162) 
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where  = 1.  For  j ^ 0.  TV.  k = 0,  we  have  zq  = 1 and.  by  L'Hospital’s  rule, 


(163)  g\  (z0)  = 


CjlV2 


-Z-  lim 
1 - 


sin  (TV  cos  1 2)  \ A/-2 


(1- 


,2^1/2 


1 - 


co(-l)J 
cj  (!  - zj)' 


where  cq  = 2.  For  j ^ 0,  AT,  = TV,  we  have  2 at  = — 1 and.  by  L'Hospital's  rule, 


(164)  g’Ax  N) 


(~1)J+1 

cjN2 


N 


1 + 


lim 


sin  (TV cos  1 z)  ^ TV2(  — 1)N 


- ,2^1/2 


(I-22) 


1 Zj 


CN(-iy+N 

Cj  (zn  - Zj) 


where  cyv  = 2.  For  j = 0,  k ^ 0.  TV, 
(165) 


5o(2fc)  caN2  ^(1  + z^TN(zk)}  Co(2fc_!) 


where  = 1,  cq  = 2.  For  j — 0,  k = 0 we  start  with 


(166) 
so  that 
(167) 


9o(z)  — 7T^2  [(1  + z)T'n(z)\ 


9o(z)  — 7^2  1^n(z)  + (1  + 2)7jv(2)]  • 


Since  g'0  (20)  = lim^-^i  g'Q(z)  we  need  to  find  2^(1)  and  T^(l).  From  the  construction  of 
T'n(z)  and  L’Hospital’s  rule, 


(168) 

Also 


2^(1)  = —TV  lim 


z — y 1 


sin  (TV  cos  1 2) 
(1-22)1/2 


= TV2 


T'M  = —TV  lim 

Z — >■  1 


TV  (l  — 22) 1//_  cos  (TV  cos  1 2)  + 2 sin  (TV  cos  1 2) 
(1  -22)3/2 


(169) 

Therefore 

(170) 

For  j yk  0,  TV,  we  use 

(171) 


TV(l  — TV2)  /sin  (TV cos  1 2)  \ TV4  — TV2 
— - hm  1 — 1 — 


z — y 1 


9o  ( zo ) = So(!)  = 


(1-22)1/2 

2TV2  + 1 


^ (zj) 


T'n  (zj)  = 


( — 1)J+1TV2 
I-22  ’ 

3(-1)j'  + 1TV2Zj 


c,  = 1,  and  L’Hospital's  rule  to  show 


9j  (Zj)  = — — — lim 


(172) 


TV2 

(-DJ+1 

2N2 


-2 zTfa)  + (1-^2)T"(2)  (1-22)T^(2) 


(2  - Zj)  ( Z-Zj ) 

[-42jr"(2J)+(l-2J2)T"'(2J)]  =- 


2(1  -ZjY 


Finally,  for  j = TV,  k = N,  cn  = 2, 
(173) 


Viv  for)  = .^2  J_imi  [-^(2)  + (1  - z)T^(2)] 
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By  L’Hospital’s  rule 


, sin  (TV  cos  1 2) 

T'n  -1  = -TV  Hm  i —L 

z~>~1  (1-22)V2 


(174) 

Also,  by  L’Hospital’s  rule, 


T'n(-I)  = Hm 

Z— > — 1 


-TV2(-1) 


— TV2  (l  — z2)1^2  cos  (TV  cos  1 2)  — TVz  sin  (TV  cos  1 z) 


N 


N3  - TV 


lim 


(175) 
Therefore 

(176) 


3 z->-i 

TV4  - TV2 


(1  - z2)3/2 
sin  (TV  cos-1  2) 


2\!/2 


:-i) 


AT 


(l-z2) 


2 TV2  + 1 


Sat  (zn)  — 7.  — ~~9o  (^0)  • 


16.  Appendix  2:  Bounds  and  Lipschitz  Condition  for  R(z,(3) 

In  this  section  we  give  a proof  of  Lemma  4.1.  A lengthy,  but  direct,  calculation  shows 


R(z,P)  = 


(177) 


From 

(178) 

we  have 

(179) 
Similarly 

(180) 

Also,  from 


(181) 
we  have 
(182) 


I / \ CcJ  \ 

Xi  \x  + s 2,  xqj  + s (x^  - xcj)  4-  s—z^  - Xx  {x,xc>) 
uj  u> 


— zds 

UJ 


- 1 


+ 


+ 


X2  ( X + s—z,  Xi  + S ( X u,  - Xi)  + S — Z u ) - X2  (x,  Xi) 
UJ  u 

X2  X + S 2,  Xjj  + s {Xu  - Xi)  + s — zw  - X2  (x,  Xi) 

(jJ  UJ 


V 1 
z^as 
U 


(iu,  - Xjj)ds 


+ ( - - 1 Xi  (x,  X&)  2 + 1 X2  (x,  Xi) 


UJ 


+ \x2  (x, Xi)  (x^  - Xi)  4-  /3X2  (x, Xi)  xc 

UJ 

+ -Xi  (x,Xi)  {zw  - Zi) . 


= f* 


(t  — uj  — s/3 ) {—(3)ds 


\Xuj  Xi  — 


| | hi  I [3 1 | Z | . 


(iu  — Xi)  + (3x i = — (3  / Xi  (t  — a>  — s/3)  — Xi  (t  — a>)  ds 


./o 


1 />! 


= 01 


0 ./o 


x(t  — lj  — us/3)  s du  ds 


(32  I- 


xw  - Xi)  + /?Xi|  < — |x| 
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Using  (177)  through  (182),  along  with  (3)  and  (4),  we  have 
(183)  \R(z,(3)\<K(z,(3), 


where 


(184) 


K(z,(3) 


2/C  I-  \z\2  + 2^-y-  (|co|  + \/3\) 

I ^ |co| 

+ — f /C  | x | + B | x j 'j  + B 


To  establish  the  Lipschitz  condition  we  start  with  the  inequality 


| dX  (x  + aii,x&  + ai2;&n,  &12)  - dX  (x  + a2 1,  + a22;  &21,  ^22)! 

(185)  < /C  (|6n|  + |5i2|)  (|an  — Q-21 1 + |o  12  — a22|) 

+/C  (|ai2|  + |a22|)  (|6n  — fe2i|  + |5i2  — &22|) 

+B  (|6n  — fe2i|  + |i>i2  — 522|) 


We  need  to  define  some  functions  that  will  help  simplify  the  relations  somewhat.  Let 


(186) 


7 = s/3  + (1  — s)0 
q = sz+(l-s)z 
9 = sz  + {\-s)z 


for  0 < s < 1.  and  define 


(187) 


0i(9,7) 

02(9,7) 

0i(9,7) 

02(9,7) 


u) 

x H q 

CO  + 7 


— xCj+-i  + 


co  + 7 


9u>+7 


(Co  + ^y 


00 


to+7  7 9^+7 

co  + 7 


CO 


(u)  + 7) 


r 9u>h — > 


Since  we  have  earlier  chosen  /3,  (3  so  that 


(188) 

it  is  easy  to  see  that 
(189) 


^ > f 

d) 

^ + /3  > ^ 


Co  + 7 
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From  (186)  we  have  the  following  integrals 


J \q\2ds  < -(|z|  + |5|) 

J h\ds  < i (|/3|  + j/3 

i 


\q\da  < ~{\z\  + \z\ 


(190) 


P\ ) (M  + l2l) 


j h\\q\ds  < 3 

I \q\\q\ds  < 7 (1*1  + |*[)  (M  + \z\) 

j IqIItMs  < \ (|i|  + \z\)  (|/5|  + |/?|) 

f1  1 

j \<i\ds  < 2 (1^1 + l5D 


Define  the  function 


:i91) 


u J UJ 

F{z,0)  = X x + rz,xo+i3  -1 — ~jZx+t3 

X ~r  p X p 


Taking  partial  derivatives  of  (191), 


UJ  ( UJ  UJ 

diF(z,(3;  y)  = — —dX(x+ — — z,x&+p  + — — zc,+p-,y,yc,+p 

l X + 0 \ X + 0 X + p 

(UJ  u) 

£ + — —32,  xcj+0  H Ftj+p; 

X + p (X  + p 


(192) 


x 


X 


£ui+/3  T 

X + p 


X 


(U  + 0) 


n\2 


•7  -uj+/3 


(u  + PY 

d\F(0,  0;  y)-  = dX  (x , x&;  y,  yU) 
d2F(0,  0;  Tj)  = r)dX  (x,  x 0,  —x^ 

From  the  definition  of  R(z,0)  and  (191)  we  have 

(193)  R(x,  0)  — R (z,  p)  = F(z,  0)  — F (z,  p)  - d2F  (o,  0;  P - p)  - d,F  (0.  0;  * - 
From  the  definition  of  7 and  q in  (186)  we  define  the  derivative  with  respect  to  s as 

(194)  dsF(q,  7;  ds)  = \d\F  ( q , 7;  2 — z)  + d2F  (q,  7;  0 — 0j  j ds. 

By  the  Fundamental  Theorem  of  Calculus 

(195)  J dsF(q,  7;  ds)  = F(z,  P)  — F [z,  0^J 
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We  can  the  write,  using  (187)  and  (192) 

R{z,0)  — r(z,i 


[d\F  (q,  7;  2 — z)  — d\F  (0.  0 ; 2 — z)\  ds 


(196) 


-i 

+ J \d2F  (9,7;/?  - (3)  - d2F  (o,0;/?  - /?)j  ds 

= / [dX  {ipi(q,'y),Tj>2{q,'v)',il’i{z,K)  - 0i(2,7),02(2,7)  - 02(2,7)) 
do 

-dX  (-01  (0, 0),  02(O, 0);  -01  (2, 0)  - 0i(5, 0),^2(z,0)  - 02(i,O))]  ds 


dX 


(■01(9,7), ^2(9, 7);  (/3  - /?)  0i(9,7),  (/?  - /?)  02(9,7)) 

—dX  (1/7(0, 0),02(O, 0);  (/3  - /?)  0i(0, 0),  (/?  - 0)  02(0. 0))]  ds 

From  (187)  we  note  that  1/7  (0,0)  = x and  02( 0,0)  = x&. 

Then,  using  (185)  through  (196)  it  is  possible  to  show  with  some  effort  that 


Xi  ( 2 


(2,/?,  2,/?)  = 


= 8/C  (I2I  + 151)  + 2/C  x + 


5 


M 


+ 0 


7^2  (2,^,2,^)  — — (—  + 4)  (|2|  + |5 


(197) 


+ ^ 2/C  | X | + 5 ^ 1 + yry  J J ( 1 2 1 + | 2 | 

2B  I -IN 
+Td  + 

\Lu\ 

i6/C  ..  . ,,,, 

+tti  (N  + M)OzI  + M) 

O (U 


B x 


+- 


+ 


H) 


17.  Appendix  A. 3:  Bounds  and  Lipschitz  Conditions  for  S{g) 

Let  g € AC  and  let  r = 5.  Then  from  Lemma  5.1  and  the  selection  of  0 so  that  F + 0(g)  > 
j,  we  have 


(198) 

and 


(199) 


u> 


^ + 0(g) 


< 2 


\S(g)\  = \Ho(z(g),0(g))\ 

2 

1,) 

< 2/C 


+2 


F + 0(g) 

l/?(g)l|2(g)| 

|0>  + /?(g)l 


l*(g)l2 

m + \P(g)\) 


+ 


0(gf 


I . i2  1. | 

/C  x + B x 


+B 


u> 


F + 0(g) 


\0(g)\\z(g)\ 
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If  we  combine  (40),  (198),  and  (199)  we  have 

\S{g)\  < J 320^  + l(,V°.Al  (pi  + 2X0S) 


(200) 


Set 


I6A0A1 
|u)| 

+2A02  ( K |f  j 2 + B if  I ) + SBX0X2  ^ 


Ei(6)  = \ 32/CAi2  + (p|  + 2X06) 


(201) 

and  let  Fi  be  a positive  constant  such  that 

(202) 


I6A0 A] 

Cj\ 

+2A02  ( K If  T + B If  I ) + 8£A0A2  U 


F\  > 32/CAi2  + 1(,A('Al  (M  + 2A0<5) 


+2A0  ^/C  j f I +B|f|l  +8£A0A2. 

Now  let  g,  5 G J\f  and  again  set  r = S.  Then,  from  (40),  (46),  and  (197)  and  .choosing 
M < <5-  we  have,  with  some  algebra, 


\S(g)-S(g)\  < 


Ai<  8/C  (|z(ff)|  + \z  (5)!) 


B 


+ ( 2/C  I f | -4  — ) {\0{g)\  + |/3  (5)!) 


+^0^  — ^ 7-7  + 4 ) .(|2(5)|  + 1 2 (5)!) 


+ (^2/C  |f  | + 13  ^ 1 -+  — j j (|2(5)|  + \z  (5)!) 
+7 — 7 (|2(5)l  + \z  (5)1) 

I UJ  I 

+ 7^(M5)l  + |2(5)l)(|i(5)l  + |i(5)l) 

O \Ld 


B f 

+-41  (|/%)l  + 1/3  (5)1) 


(203) 


< 


Ai<(  32/CAi  + 4A0  ( 2/C  f + — 


\9~9\ 


B 


+Aq 


16A02  / 16 


— +4  I <5  + 4Ai  ( 2/C|f| +/3  ( 1 + — 


8BX2  256/C  AiA9r 
^ + o-r-7--d 

IlJ  I o \uj\ 
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Finally,  we  set 

E2(6)  = 

(204) 


A}  ^ 32/CAi  + 4Aq  ( 2/C  x j + ■ 


+Aq 


16A02  ( 16 


M 


— + 4 )5  + 4Ai  2/C \x  + B 1 + — 


8BA2  256/CAxA2  p 
H — rrr-  + — - , . . "6 


3 M 


+2B  x Af 


and  let  !■'■>  be  a positive  constant  such  that 

F2  > Ai  32/CAi  + 4A0  f 2/C  If  I + ~ 


(205) 


| 16V  /_16 


+Aq 


8B\2  256/CAiA2  r 
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18.  Appendix  3:  Main  Matlab  Script 

This  section  includes  the  main  script  and  supporting  functions,  except  for  “cheb.m”, 
which  is  available  in  Trefethen  [31].  These  scripts  are  included  as  is.  They  are  not  necessarily 
the  most  efficient  and  are  specifically  oriented  towards  the  Van  der  Pol  equation  example  in 
Section  12.  A user  will  have  to  modify  the  scripts  for  their  particular  problem, 
global  m N CS  MO  M2  VO  VI  V2  T lambda  DM 
global  a_bar 
global  startt  endt 
global  D_hat 
global  A.. hat 
global  NC 
global  zj 
global  m 
global  piinvN 

y,****************************************************"^ 

%User  input \\ 

%*************#*****:i:>l=*=i'***=t<=i=*********=i=*>i'*>l'**********\\ 

cal_B  = input ( ’Bound  on  derivatives  for  right  hand  side  of  DDE.  cal_B  = ’); 
cal_K  = input ( ’Lipschitz  condition  on  right  hand  side  derivatives  of  DDE.  cal_K  = 
lambda  = input(’Van  der  Pol  Equation  parameter  lambda,  lambda  = ’); 
m = input (’ Approximate  Solution  Harmonics,  m = ’); 

N = input ( ’ 1/2  number  of  integration  points  for  init . cond.  function.  N = ’); 

NC  = inputf’Enter  NC  for  NC  + 1 Chebyshev  points  for  collocation.  NC  = ’); 

P = input (’Enter  number  of  points  for  Trapezoidal  integration  on  [-omega, 0] . P = ’ 
0 = input(’Enter  number  integration  points  on  [0,2pi] , P/0  = 1/6,  e.g.  250/1500.  0 


^Computing  the  initial  approximation  fnnctionW 

%5^5|e^^:^S5je5ie3jc3jc5^^3^^^e^e3(c3^j|:^5f£^3^3is3^^:^^c^^e5i:^:^3(s3ic^3f:3ic^:^3|c^3is5tc3f:*3|esjc3^^5|c^:3f:3|c3i:W 
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’/initialize  arrays 
T = zeros (2*N, 1) ; 
i = zeros (2*N , 1) ; 

VO  = zeros (2*m-l , 1) ; 

VI  = zeros(2*m, 1) ; 

V2  = zeros(2*m, 1) ; 

MO  = zeros(2*N,2*m-l) ; 
M2  = zeros(2*N,2*m-l) ; 
A = zeros (2*m , 1) ; 

AO  = zeros(2*m, 1) ; 

DM  = zeros(2*N,2*m) ; 

CS  = zeros(2*m,2*N) ; 


’/Open  I/O  file 

fid  = f open( ’Est_Periodic_Sol_err_Output . txt ’ , ’w’ ) ; 


’/Projection  integration  steps 
i = (1 : 2*N) ’ ; 

T = (pi/(2*N))*(2*i-l) ; 


’/Set  up  fixed  arrays 
for  n = 1 :m 


CS(2*n-l,:)  = cos ( (2*n-l)*T) ’ ; 
CS (2*n, : ) = sin( (2*n-l) *T) ’ ; 


end 


MO ( : , 1)  = cos (T) ; 
M2 ( : ,1)  = -cos (T) ; 


n = 

1 :m-l 

MO  ( 

,2*n)  = 

MO  ( 

,2*n+l) 

M2  ( 

,2*n)  = 

M2( 

,2*n+l) 

cos ( (2*n+l) *T) ; 

= sin( (2*n+l)*T) ; 

(- (2*n+l) \~2) *cos ( (2*n+l) *T) ; 

= (- (2*n+l) \"2) *sin( (2*n+l) *T) ; 


end 


for  n = 1 :m 

DM(:,2*n-l)  = cos ( (2*n-l)*T) ; 
DM(:,2*n)  = sin( (2*n-l) *T) ; 

end 


’/Initialize  AO 
A ( 1 ) = 1.0; 

A (2)  = 2.0; 

a_bar  = f solve( ’Numerical.Galerkin’ , A) 

fprintf (fid, ’Frequency  and  Approximate  Solution  Coefficients  \n$’ 
nabar  = length (a_bar) ; 
for  io  = 1: nabar 

f printf (f id, ’ %15 . 8e  \n ’ , a_bar (io) ) ; 

end 


’/calculate  the  residual 
r = Van_der_Pol (a_bar) ; 
norm\_r  = max(abs(r)) 

f printf  (fid,  ’ \n\nResidual  error  r = °/15.8e  \n’,norm_r); 
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Lgt\_T  = length (T); 
xt  = zeros(Lgt_T, 1) ; 
xdt  = zeros (2,Lgt_T) ; 
xdt\_temp  = zeros (2,1); 
for  i = 1:2*N 

xt(i,l)  = Galerkin_series (T(i) , a_bar ,m) ; 
xdt_temp  = Derivative_series (T(i) ,a_bar ,m) ; 
xdt(l,i)  = xdt_temp(l , 1) ; 
xdt(2,i)  = xdt_temp(2, 1) ; 

end 

norm_abs_xt  = max (abs (xt ( : , 1) ) ) 
norm_abs_dxt  = max (abs (xdt (1 , : ) ) ) 
norm_abs_ddxt  = max(abs (xdt (2, : ) ) ) 

fprintf (f id, ’\n\nMax  absolute  value  x(t)  on  [0,2pi]  = %15.8e\n’ ,norm_abs_xt) ; 
fprintf (f id , ’ \n\nMax  absolute  value  dx/dt(t)  on  [0,2pi]  = 7,15 .8e\n’ ,norm_abs_dxt) ; 
f printf  (f  id , ’ \n\nMax  absolute  value  d2xdt2(t)  on  [0,2pi]  = °/015 . 8e\n’  ,norm_abs_ddxt) 

“/.Plot  residual  error 

figure ; 
plot (T,r)  ; 

title( ’Residual  Error  as  a Function  of  Time  over  [0,  2\pi] ’ ) ; 

xlabel ( ’Time ’ ) ; 

ylabel (’ Residual  Error’); 

disp( ’press  enter  to  continue’) 

pause; 

^^^^^^^^^^5je^^^^^3j«3^^:5(c3(e5t:3je^:^e^:^:^^5^^:^:jjes|c4:5ie5)c5f:5fe^£5f:5f:jjc3f:3t:5(c5(c3f:5f:^:^c^:^:5|e3jc3(e3|:^e5t:5f:3(£^e 

"/Phase  plot 

nt  = 10000; 

t = zeros(19999, 1) ; 

ii  = (1 :nt) ’ ; 

t = (ii-l)*2*pi/nt ; 

x_hat  = Phase_series (t , a_bar ,m) ; 

figure; 

plot (x_hat ( : , 1) ,x_hat ( : ,2) ) ; 

title (’Phase  Plot  of  Approximate  Solution’); 

xlabel ( ’ X ’ ) ; 

ylabel (’ dX/dt ’) ; 

axis  equal; 

disp( ’press  enter  to  continue’) 
pause ; 

% End  of  appoximation  section 

7. 

7o  Output  from  this  section  a_bar,  norm_r 

%*************************************************************** 

70  Begin  collocation  section 

^3|c^c3te3f:^:^c^c3f:^c^c3(e^:3f;^:^})e^e^c^c^e3f;3{c]je3|c^e^e}f;^;3(c4;3|e^c3|c^c](c^e^c3)c^c3(c^e^:^e}(e3|c3(e3(c^c3f:^e^c3f:3f:^c3f:^3f:^;^c3|c^3(c3f: 

%*************************************************************** 

1 Get  NC  Chebyshev  points,  z j , and  Create  Chebyshev  differentiation 
7.  matrix,  D.  See  Trefethen  book. 
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%*************************************************************** 
[D,zj] =cheb(NC) ; 

%*************************************************************** 
‘/»  Reset  last  row  of  D to  account  for  initial  condition  row 
% *************************************************************** 
D (NC+1 , 1 : NC)  =0.0; 

D(NC+1,NC+1)  = 1.0; 

D_hat  = zeros (2* (NC+1) ,2* (NC+1)) ; 

D_hat (1 : NC+1 , 1 :NC+1)  = D(1 :NC+1 , 1 :NC+1) ; 

D_hat (NC+2 : 2* (NC+1) ,NC+2 : 2* (NC+1) ) = D ( 1 :NC+1 , 1 :NC+1) ; 


%***************************************************************** 

y.  Set  up  fixed  arrays  for  the  Van  der  Pol  problem 

°l  Create  the  A_hat  matrix  in  the  Van  der  Pol  case.  It’s  constant. 

% Initialize  the  B array  and  a temporary  array 

y^^c3(c3(c^c^:5(c3|c^:5f:3f:3(c^c3|c^c5(c5)c^c^cjf:3|c5(c^c^:^e^e3)c^c5f:3t:^:5t:3f:^c3jc  + ^:^:3|c^c^c3f:5(e^c  + ^c^c3t:^c3(e^c^:^cj|c3t:3tc3f:^c^c^c^c3f:^c^e^:3f: 

A_hat  = zeros (2* (NC+1) , 2* (NC+1) ) ; 

A = [0  1;  -1/ (a_bar (1 , 1) ~2)  0]; 
for  i = 1:NC 

A_hat(i,i)  = A ( 1 , 1) ; 

A_hat (i ,NC+l+i)  = A(l,2); 

A_hat (NC+l+i , i)  = A(2,l) ; 

A_hat (NC+l+i ,NC+l+i)  = A(2,2); 

end 

B = zeros(2,2);  % coefficients  of  linear  delay  term  for  Van  der  Pol 

Temp  = zeros(2,2); 

‘/.Compute  the  fundamental  matrix  at  several  points  from  0 to  2pi 
‘/.for  plotting  only 

°/0get  Lagrange  weights  from  0 to  2*pi  for  Fundamental  Solution 

disp (’ Plotting  Fundamental  matrix  from  0 to  2pi ’ ) ; 

‘/.Weights  = colloc (0 , 2*pi)  ; 
a = 0; 
b = 2*pi ; 

‘/.debug  * ***  *********  ******  *******************  ********** 

Q = fix( (b-a)/a_bar(l , 1))+1 ; 

B = zeros(2,2);  tj  = zeros (NC+1 , 1) ; g = zeros(2,l); 

W0  = zeros(2* (NC+1) , 1) ; Weights  = zeros (2* (NC+1) , 2 , Q) ; 

B_hat  = zeros (2* (NC+1) , 2* (NC+1) , Q) ; 

‘/.solve  for  weights  on  first  step  interval 
omega  = a_bar(l,l); 

M = ( (D_hat  - (omega/2) *A_hat) “ (-1) ) ; 

W0 (NC+1 , 1 ) = 1;  W0 (2* (NC+1 ) , 1)  = 0; 

Weights (: ,1,1)  = M*W0; 

W0 (NC+1 , 1 ) = 0;  W0 (2* (NC+1 ) , 1)  = 1; 

Weights ( : ,2,1)  = M*W0; 

‘/.other  intervals 
for  i = 2:Q 

B_hat (NC+1 , 1 , i)  = 2/omega;  B_hat (2* (NC+1) ,NC+2 , i)  = 2/omega; 
tj  = (omega/2)*zj  + (l/2)*(2*a  + (2*i-l) *omega) ; 
for  k = 1 : length(t j )-l 

g(l,l,l)  = a_bar (2 , 1 , 1) *cos (t j (k) ) ; 
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g(2 , 1 , 1)  = -a_bar (2, 1 , l)*sin(t j (k) ) ; 
for  n = 2:m 

g(l,l,l)  = g(l,l,l)  + a_bar(n*2, l)*cos(n*tj (k) ) + a_bar (n*2-l , l)*sin(n*t j (k) ) ; 
g(2,l,l)  = g(2,l,l)  - n*a_bar (n*2 , 1) *sin(n*t j (k) ) + n*a_bar (n*2-l , 1) *cos (n*t j (k) ) ; 

end 

B(2, 1,1)  = (-2*lambda/ a_bar (1 ,l,l))*g(l,l,l) . *g(2, 1 , 1) ; 

B (2 , 2 , 1)  = (lambda/a_bar(l,l))*(l-g(l,l,l) .~2) ; 

°/0f ill  up  the  B_hat_i  matrix 
B_hat (k,k, i)  = B(l,l,l) ; 

B_hat(k, (NC+l)+k,i)  = B(l,2,l); 

B_hat ( (NC+l)+k,k, i)  = B(2,l,l); 

B_hat ( (NC+l)+k, (NC+l)+k,i)  = B(2,2,l) ; 

end 

Weights ( : , 1 , i)  = M* (omega/2) *B_hat (1 : 2* (NC+1) ,1:2* (NC+1) , i) ^Weights ( : , 1 , i-1) ; 

Weights ( : ,2,i)  = M* (omega/2) *B_hat (1 : 2* (NC+1) , 1 : 2* (NC+1) , i) *Weights ( : ,2, i-1) ; 

end 

%end  debug  ******************************************** 

[max„row,  max_col,  max_plane]  = size (Weights) ; 

'/interpolate  the  fundamental  solution 

pp  = 1000; 

tt  = zeros (pp+1 , 1) ; 
zz  = zeros (pp+1 , 1) ; 
yy  = zeros(NC+l , 1,1); 
del  = 2*pi/pp ; 
for  ii  = l:pp+l 

%get  a time  value  betwee  0 and  2*pi 
tt(ii)  = (ii-l)*del; 

Z = vdp_interp(0 ,tt (ii) .Weights) ; 

yzll (ii , 1)  = Z(l,l) ; 

yz21  (ii , 1)  = Z(2, 1) ; 

yzl2 (ii , 1)  = Z ( 1 , 2)  ; 

yz22  (ii , 1)  = Z(2,2) ; 

end 

f igure ; 

subplot (2 ,2 , 1)  ; 

plot (tt (1 :pp+l) , yzll (1 : pp+ 1 , 1) ) ; 

title(’Z(l,  1).’ ) ; 

xlabel ( ’ tt ’ ) ; 

ylabel ( ’ zll ’ ) ; 

subplot (2, 2, 2) ; 

plot (tt ( 1 :pp+l) ,yz21 (1 : pp+1 , 1) ) ; 

title ( ’ Z (2, 1) ’) ; 

xlabel ( ’ tt ’ ) ; 

ylabel ( ’ z21 ’ ) ; 

subplot (2, 2, 3) ; 

plot (tt (1 :pp+l) , yzl2 ( 1 :pp+l , 1) ) ; 

title(,Z(l,2) ’) ; 

xlabel ( ’ tt ; ) ; 

ylabel ( ' zl2 ’ ) ; 

subplot(2,2,4) ; 

plot (tt (1 :pp+l) , yz22 (1 : pp+1 , 1) ) ; 
title( ’ Z (2,2) ’ ) ; 
xlabel ( J tt ’ ) ; 
ylabel ( ’ z22 ’ ) ; 
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disp( ’press  enter  to  continue’) 
pause ; 

l clear  out  the  arrays  not  needed 
clear  pp  tt  yy  zz  Z Weights 

•/.  Computing  the  monodromy  matrix.  This  is  a 2*  (P+1)  x 2*  (P+1)  matrix 
% Eigenvalues  by  Collocation.  Integration  by  Trapezoidal  rule. 

^sfcafcsfesfes^alcsfesfcsfcsfcafcsfcsfcsfcsfc^fcsfcafcsfcs^sfcsfcsfcsfesfcsfesfestcsfcsfcsfcsfcafcafcsfcsfcstcafcsfcsfcsfcstcafcafcsfcsfestcsfeafcsicsfcsfcsfcsfcsfcstcsfc^fcsfcsfcsfcstc  afc  sfc  sfc  afc  sfc  3fc 

omega  = a_bar(l,l); 

delta  = omega/P;  °/0  P is  user  selected  number  of  Trapezoidal  points 
iii  = 1 : P+1 ; 

s = -omega  + (iii-l)*delta;  ^integration  from  -omega  to  0 for  monodromy  matrix 
U = zeros  (2*  (P+1)  , 2*  (P+1) ) ; °/0  the  2 is  for  the  2x2  blocks  for  VdP  equat . 

Z = zeros(2,2);  'predefine  Z again 
endt  = 2*pi; 

B ( 1 , 1)  = 0.0; 

B ( 1 , 2)  = 0.0; 
deltal  = delta; 

1 Load  by  column 
for  jj  = 1:P+1 

disp  (sprintf  (’ Creating  Block  column  °/0d  for  U\n’,jj)) 
startt  = s(jj)  + omega; 

%get  the  weights  for  the  number  of  step  intervals  between 
°/„startt  = s(jj)+omega  and  endt  = 2*pi 
Weights_eig  = colloc (startt , endt) ; 

% Set  up  B matrix  for  integration  step  using  trapezoidal  rule 
x_hat  = Vdp_series (s ( j j ) ,a_bar ,m) ; 

B(2,l)  = (-2*lambda/omega) *x_hat (1) *x_hat (2) ; 

B(2,2)  = ( lambda/ omega) * (1 . 0 - x_hat(l)"2); 

% form  row  blocks  for  column  j 
for  ii  = 1:P+1 

tt(ii)  = s(ii)  + 2*pi; 

Z = vdp_interp( startt ,tt (ii) , Weights_eig) ; 
if  ((ii  ==  1)  I (ii  ==  P+1)) 
deltal  = delta/2; 

end 

Temp  = deltal*Z*B; 

% Fill  in  the  i row  block  of  the  current  column 
U(2*ii-l,2*jj-l)  = Temp(l , 1) ; 

U(2*ii-l,2*j j)  = Temp(l ,2) ; 

U(2*ii  ,2*j  j-1)  = Temp(2d); 

U (2*ii  ,2* j j ) = Temp(2,2) ; 

end 

clear  weights_eig  '(need  to  clear  since  startt  changes 

end 

%add  the  initial  condition  blocks  to  the  last  column 

disp (’Adding  the  initial  condition  blocks  to  the  last  column’) 

jj  = P+i; 

startt  = 0; 

Weights_eig  = colloc (startt , 2*pi) ; 
for  ii  = 1: P+1 

tt(ii)  = s(ii)  + 2*pi ; 
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Z = vdp_interp (startt ,tt (ii) , Weights_eig) ; 
U(2*ii-l,2*(P+l)-l)  = U(2*ii-1 ,2* (P+1) -1)  + Z (1 , 1) ; 
U(2*ii,2*(P+l)-l)  = U(2*ii,2*(P+l)-l)  + Z(2,l) ; 
U(2*ii~l,2*(P+l))  = U(2*ii-l,2*(P+l))  + Z ( 1 , 2) ; 
U(2*ii , 2* (P+1) ) = U(2*ii,2*(P+l))  + Z(2,2); 


end 

eig(U) 

disp(’ Operator  matrix  U is  filled,  now  computing  eigenvalues’) 

[V.Diag]  = eig(U) ; 

[row, col]  = size(Diag); 
abs_diag  = zeros (row, 1) ; 

disp( ’Eigenvalues  of  U by  colloc  and  max  absolute  value’) 
fprintf (f id, ’ \n\nEigenvalues  of  U\n’); 
for  i = 1 :row 

eigenCi , 1)  = Diag(i,i); 
re  = real (eigenCi , 1) ) ; 
im  = imag( eigenCi , 1) ) ; 

fprintf (fid, ’ °/„15 . 8e  + i °/,15 . 8e\n’ ,re , im) ; 

end 

for  i = l:row 

abs_diag(i , 1)  = abs (Diag(i , i) ) ; 

end 

max_abs_diag  = max (abs_diag) 

fprintf (fid, ’ \n\nMaximum  Absolute  Value  of  Eigenvalue  for  U\n’) 
fprintf (fid, ’%15 . 8e\n’ ,max_abs_diag  ) ; 

%plot  the  first  20  eigenvalues  of  U 

figure 

hold  on; 

ang  = 0 : pi/100 : 2*pi ; 
plot (sin (ang) , cos (ang) , ’b- ’ ) 
plot (eigen, ’rx’ ) 
hold  off; 

title (’ Eigenvalues  for  U’); 
xlabeM’Real  Part  of  Eigenvalue’); 
ylabel (’ Imaginary  Part  of  Eigenvalue’); 
axis  equal; 

°/.disp( ’press  enter  to  continue’) 

"/.pause ; 

%********************************************************************* 
"/.Computing  solution  of  the  adjoint  associated  with  characteristic 
"/.multiplier  near  the  unit  circle 

%********************************************************************* 
"/.First  fill  the  matrix  U_tilde.  According  to  Halanay  this  should  have 
"/.the  same  eigenvalues  as  U.  Trapezoidal  integration  from  -omega  to  0. 

y********************************************************************* 
U_tilde  = zeros (2* (P+1) ,2* (P+1) ) ; 
delta  = omega/P; 
iii  = 1 : P+1 ; 

s = -omega  + (iii-1) *delta; 
deltal  = delta; 

B (1 , 1)  = 0.0; 

B(1 ,2)  = 0.0; 

"/.load  by  columns 


38 


DAVID  E.  GILSINN 


for  jj  = 1 : P+1 

disp (sprintf (’ Creating  Block  column  %d  for  U_tilde\n’ , j j ) ) 
startt  = s(jj)  + omega; 

Weights_adj_eig  = colloc (startt , 2*pi) ; 

"/.Load  row  blocks  for  colum  jj 
for  i i = 1 : P+ 1 

%B  Matrix  for  Van  der  Pol 

x_hat  = Vdp_series (s (ii) , a_bar ,m) ; 

B(2,l)  = (-2*lambda/omega) *x_hat (l)*x_hat (2) ; 

B(2,2)  = ( lambda/ omega) * (1 . 0-x_hat (1) ~2) ; 

"/.Interpolate  at  row  ii  Z(s(ii)+2pi  ,s(j  j)+omega) 

Z = vdp_interp(startt,s(ii)+2*pi , Weights_adj _eig) ; 
if  ((ii  ==  1) I (ii  ==  P+1)) 
deltal  = delta/2; 

end 

Temp  = deltal*B*Z; 

70Fill  in  the  row  blocks  ii  of  the  current  column  jj 
U_tilde (2*ii-l , 2* j j -1)  = Temp(l,l); 

U_tilde (2*ii-l , 2* j j ) = Temp(l,2); 

U_tilde(2*ii ,2*j j-1)  = Temp(2,l); 

U_tilde (2*ii , 2* j j ) = Temp(2,2); 

end 

clear  Weights_adj_eig 

end 

"/.add  initial  condition  block  to  first  row 

disp( ’Adding  the  initial  condition  blocks  to  the  first  row’) 
endt  = 2*pi ; 
for  jj  = 1 : P+1 

startt  = s(jj)  + omega; 

Weights_adj_eig  = colloc (startt , endt) ; 

"/.Interpolate  Z(2pi , s(  j j )+omega) 

Z = vdp_interp (startt , endt ,Weights_adj _eig) ; 

U_tilde(l ,2*j j-1)  = U_tilde (1 ,2*j j-1)  + Z(l,l); 

U_tilde(2,2*j j-1)  = U.tilde (2 ,2* j j-1)  + Z(2,l); 

U_tilde ( 1 , 2* j j ) = U_tilde(l ,2*j  j ) + Z(l,2); 

U_tilde (2 , 2* j j ) = U;_tilde(2,2*j  j)  + Z(2,2); 
clear  Weights_adj _eig 

end 

eig(U_tilde) 

disp (’ Operator  matrix  U_tilde  is  now  filled,’); 
disp (’ computing  eigenvalues  and  eigenvectors  of  transpose’); 
disp ( ’Transposing  U_tilde  is  necessary  in  order  for  Matlab’); 
disp(’to  get  column  eigenvectors’); 

[PhiT.Dadj]  = eig(U_tilde’ ) ; 

[row, col]  = size(Dadj); 
abs_diag  = zeros (row , 1) ; 

disp ( ’Eigenvalues  by  Colloc  of  U-tilde  transpose  and  max  abs  val. 
f printf (f id , ’ \n\nEigenvalues  of  U-tilde\n’); 
for  i = l:row 

eadj(i,l)  = Dadj(i,i); 
re  = real (eadj (i , 1) ) ; 
im  = imag(eadj (i , 1) ) ; 

f printf  (fid,  ’ ”/„15 . 8e  + i 7,15 . 8e\n’  ,re  , im)  ; 


end 
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for  i = 1 :row 

abs_diag(i , 1)  = abs (Dadj (i , i) ) ; 

end 

max_abs_diag  = max (abs_diag) ; 

fprintf (fid, ’ \n\nMaximum  Absolute  Value  of  Eigenvalue  for  U-tilde\n’) 
fprintf (fid, ’“/.15 . 8e\n’ ,max_abs_diag  ) ; 

I_max_eig  = f ind(abs_diag  ==  max_abs_diag) 

figure; 

hold  on; 

ang  = 0:pi/100:2*pi; 

plot (sin (ang) , cos (ang) , ’b-  ’ ) ; 

plot (eadj , ’rx’) ; 

title ( ’Eigenvalues  of  U-tilde’); 
xlabel ( ’Real  Part  of  Eigenvalue’); 
ylabel (’ Imaginary  Part  of  Eigenvalue’); 
axis  equal; 

“/.disp( ’press  enter  to  continue’) 

“/.pause ; 

%***************************************************************** 

“/.The  first  eigenvector  in  PhiT  is  associated  with  the  first  eigenvector  of 
“/.largest  magnitude.  Transpose  it  to  make  it  a row  vector. 

pp( : , 1)  = PhiT ( : , I_max_eig(l) ) ; 

Igth  = length (pp); 

phi  = zeros (1 , Igth) ; 

phi (1, 1 : Igth)  = pp ( : , 1) ’ ; 

U_adj  = zeros (2* (P+1) ,2*(0+l)) ; 

°/.N+l  blocks 
iii  = 1 : P+1 ; 

s = -omega  + (iii-l)*delta; 

sp2pi  = s + 2*pi; 

endt  = 2*pi ; 

kk  = 1:0+1; 

t = (kk-1) *2*pi/0 ; 

for  k = 1 : 0+1 

“/.Compute  adjoint  solution  at  the  t(k) 

“/.This  will  be  the  k-th  column 

disp(sprintf  ( ’Computing  adjoint  at  y (t  (7.d) ) \n’ ,k) ) ; 

“/.Load  up  the  column  k 
startt  = t (k) ; 
if  (k  ==  1) 

startt  = 0.0; 

end 

“/.Get  weights  for  Z(2pi,t(k)) 

Weights_adj  = colloc(startt,endt) ; 
delal  =. delta; 

B(1 , 1)  = 0.0; 

B(1 ,2)  = 0.0; 

“/»  Load  row  block  for  column 
for  ii  = 1 : P+1 

%disp(sprintf  ( ’Creating  Block  column  “/.d  for  adj oint\n’ , ii) ) 
x_hat  = Vdp_series (s (ii) , a_bar ,m) ; 

B(2,l)  = (-2* lambda/omega) *x_hat (1) *x_hat (2) ; 

B(2,2)  = (lambda/ omega) * (1 . 0-x_hat (1) ~2) ; 
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'/Check  relation  between  start  and  end  time  of  fundamental  solution 
if  (startt  > sp2pi(ii)) 

Z = zeros(2,2) ; 
elseif  (startt  ==  sp2pi(ii)) 

Z = eye (2) ; 

else 

'/Interpolate  Z(s(ii)+2pi ,t(k)) 

Z = vdp_interp(startt , sp2pi (ii) , Weights_adj ) ; 

end 

if  ((ii  ==  1)  I (ii  ==  P+1)) 

deltal  = delta/2;  '/Trapezoidal  rule  end  weights 

end 

Temp  = deltal*B*Z ; 

'/  Fill  in  the  row  blocks  of  the  current  column 
U_adj (2*ii-l , 2*k-l)  = Temp(l,l); 

U_adj (2*ii-l , 2*k)  = Temp(l,2); 

U_adj (2*ii ,2*k-l)  = Temp(2,l); 

U_adj (2*ii , 2*k)  = Temp(2,2); 

end 

clear  Weights_adj 

end 

endt  = 2*pi ; 

'/add  initial  condition  block  to  first  block  row  of  column  k 
for  k = 1:0+1 

startt  = t(k); 

'/Interpolate  Z(2pi,t(k)) 
if  (k  ==  1) 

startt  = 0.0; 

end 

Weights_adj  = colloc(startt ,endt) ; 

Z = vdp_interp (startt , endt , Weights_adj ) ; 

U_adj(l,2*k-1)  = U_adj (1 ,2*k-l)  + Z ( 1 , 1) ; 

U_adj (2,2*k-l)  = U_adj (2,2*k-l)  + Z ( 2 , 1 ) ; 

U_adj (1 , 2*k)  = U_adj ( 1 , 2*k)  + Z(l,2); 

U_adj (2 , 2*k)  = U_adj(2,2*k)  + Z(2,2) ; 
clear  Weights_adj 

end 

'/Apply  row  eigenvector  to  U_adj 

°/v_0  is  1 row  by  2*  (0+1)  columns  since  phi  is  1 row  by  2*  (P+1)  columns 
'/and  U_adj  is  P+1  rows  by  0+1  columns 
°/v_0  might  be  complex 
v_0  = phi*U_adj ; 

'/normalizing  v_0 
norm_v_0  = norm(v_0,2); 
v_0  = v_0/norm_v_0 ; 

'/Computing  the  alpha  parameter 

'/  Now  compute  the  alpha  parameter  in  Stokes  theorem 
disp( 'Computing  alpha’); 

'/  First  compute  the  J (omega, x_dot)  function  at  the  mesh  points  in  [0,2+pi]. 
'/  compute  alpha.  First  comput  J(i) 

B(1 , 1)  = 0.0; 
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B ( 1 , 2)  = 0.0; 
alpha  = 0.0; 

J1  = zeros(2 ,0+1) ; 

Jl_temp  = zeros (2 , 1) ; 
w_trap  = zeros (0+1 , 1) ; 
for  k = 1:0+1 

w_trap(k,l)  = 2*pi/0; 
if  ( (k==l) I (k==0+l)) 
w_trap(k,l)  = pi/0; 

end 

'/first  set  up  B(t(k)) 

x_hat  = Vdp_series (t(k) -omega, a_bar ,m) ; 

B(2,l)  = (-2*lambda/omega)*x_hat(l)*x_hat(2) ; 

B(2,2)  = (lambda/ omega) * (1 .0-x_hat(l) ~2) ; 

'/, Multiply  by  the  derivative  at  t(k)  - omega 
x_hat_dot=delay  = Derivative_series (t (k) -omega, a_bar ,m) ; 

°/0add  derivative  at  t(k) 

x_hat_dot_t  = Derivative_series (t (k) , a_bar ,m) ; 

'/compute  J as  2 rows,  1 column 
J = x_hat_dot_t  + B*x_hat_dot_delay ; 

J_n(k,l)  = norm(J,2); 
v0(l , 1)  = v_0(l,2*k-l) ; 
v0(l ,2)  = v_0(l,2*k); 

'/form  alpha 

alpha  = alpha  + w_trap(k, 1) *v0*J ; 

end 

alpha  = 1/alpha 
norm_J  = max(J_n) 

f printf (fid, ’\n\nPrinting  estimate  of  norm  of  J on  [0,2pi]  = °/15 . 8e\n’ ,norm_J) 
f printf (fid, ’ \n\nPrinting  alpha\n ’ ) ; 
re  = real (alpha); 
im  = imag (alpha) ; 

fprintf  (f  id, ’°/15.8e  + i °/15 . 8e\n’ ,re , im)  ; 

'/This  alpha  might  be  complex  but  we  really  only  use  abs  (alpha) 
disp( ’Absolute  value  of  alpha’); 
abs _ alpha  = abs (alpha) 

fprintf (f id , ’ \n\nAbsolute  Value  of  alpha\n’); 
f printf  (f  id , ’ '/15 . 8e\n’  , abs_alpha)  ; 

°/0disp( ’Press  Enter  to  Continue’); 

'/pause ; 

'/Now  compute  the  M bound 

^ ^ 5^  s(c  3(c  sf:  sfc  3(c  sfc  3(c  5+:  3^:  sf:  sfs  3|c  ^ sfc  sf:  sfe  sf:  sf:  af:  =|c  sf:  s|c  sfc  sf;  sfc  3<c  sfc  54c  sfc  sf:  sfc  3{c  sfc  afc  a(c  af: 

'/  get  integral  discretization  steps 

ds  = omega/P; 
dt  = 2*pi/0; 

'/  Identify  points 

iP  = 1 : P+1 ; 

sP  = -omega  + (iP  -l)*ds; 
iO  = 1:0+1; 
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tO  = (iO 

"/Trapezoidal  weights  for  integration 

w = zeros (P+1) ; 
u = zeros (0+1) ; 
w(l:P+l)  = ds ; 
u ( 1 : 0+1 ) = dt; 
w ( 1 ) = ds/2 ; 
w(P+l)  = ds/2; 
u(l)  = dt/2 ; 
u(0+l)  = dt/2; 

7.  Set  up  some  fixed  matrices 
% Get  collocation  weights  between  0 and  2pi 
Weightsl  = colloc(0,2*pi) ; 

7,  Interpolate  to  form  fundamental  solution  at  2pi,  i.e.  Z(2*pi,0) 
Z = vdp_interp(0,2*pi ,Weightsl) ; 

% Form  (i  - Z(2*pi ,0) ) ~ (-1) 

12  = eye  (2) ; 

ImZ_inv  = (12  - Z) “ (-1) ; 

7,  Form  (i  - Z (2*pi ,0) ) ~ (-1) *Z(2*pi , 0) 

ImZ_inv_Z  = ImZ_inv  * Z; 

7.  Set  up  some  cell  arrays 
disp( ’ Setting  up  cell  arrays’); 

B ( 1 , 1)  = 0; 

B(1 ,2)  = 0; 
for  iPP  = 1:P+1 

Z = vdp_interp(0 , sP (iPP)+2*pi , Weightsl) ; 

Z_cell (iPP , 1)  = {Z*ImZ_inv>; 
x_hat  = Vdp_series (sP (iPP) , a_bar ,m) ; 

B(2,l)  = (-2*lambda/omega) *x_hat (1) *x_hat (2) ; 

B(2,2)  = (lambda/omega) * (1 . 0 - x_hat(l)~2); 

B_cell(iPP,l)  = {B>; 

end 

clear  Weightsl; 

7»  Initialize  the  H arrays 
IP  = eye (2* (P+1)) ; 

HI  = zeros(2* (P+1) , 2* (P+1) ) ; 

H2  = zeros(2* (P+1) , 2*(0+l)); 

H3  = zeros(2* (0+1) , 2* (P+1) ) ; 

H4  = zeros(2* (0+1) , 2*(0+l)); 

H5  = zeros(2* (0+1) , 2*(0+l)); 

% Compute  HI 
disp(’ Computing  HI’); 
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for  jl  = 1:P+1 

disp(sprintf ( ’Column  block  %d  of  %d  columns  of  Hl\n’ , j 1 ,P+1) ) ; 
V/eighf s2  = colloc(sP(jl)+omega,2*pi)  ; 

Z1  = vdp_interp(sP( jl)+omega,2*pi , Weight s2) ; 
for  il  = 1 :P+1 

Tempi  = Z_cell{il , 1}*Z1 ; 

Z2  = vdp_interp (sP ( j 1) +omega , sP(il)+2*pi , Weights2) ; 

Temp 2 = (Tempi  +Z2) *B_cell{ j 1 , 1} ; 

Temp 2 = w(jl)*Temp2;  ’/.Trapezoidal  rule  weights 
Hl(2*il-l,2*jl-l)  = Temp2(l , 1) ; 

HI (2*il-l ,2* j 1)  '=  Temp2 (1 , 2) ; 

HI (2*il , 2* j 1-1)  = Temp2 (2,1); 

HI (2*il , 2* j 1)  = Temp2(2 , 2) ; 

end 

clear  Weights?. ; 

end 

7.  Compute  H2 

disp( ’Computing  H2’); 

for  k2  = 1:0+1 

disp(sprintf  ( ’Column  block  °/.d  of  7,d  columns  of  H2\n’  ,k2,0+l) ) ; 
WeightsS  = colloc (t0(k2) , 2*pi) ; 

Z1  = vdp_interp(t0(k2)  , 2*pi , Weights3)  ; 
for  i2  = 1 :P+1 

if  (tO (k2)  > sP(i2) +2*pi) 

Z2  = zeros(2,2) ; 
elseif  (t0(k2)  ==  sP(i2)+2*pi) 

Z2  = eye (2) ; 

else 

Z2  = vdp_interp (tO (k2) ,sP(i2)+2*pi ,Weights3) ; 

end 

Temp 3 = Z_cell{i2 , 1}*Z1  + Z2; 

Temp 3 = u(k2)*Temp3; 

H2(2*i2-1 ,2*k2-l)  = Temp3(l,l); 

H2(2*i2-1 ,2*k2)  = Temp3(l,2); 

H2(2*i2,2*k2-1>  = Temp3(2,l) ; 

H2(2*i2 , 2*k2)  = Temp3(2,2); 

end 

clear  WeightsS; 

end 

Weightsl  = colloc (0 , 2*pi) ; 

% Compute  H3 

disp( ’Computing  H3 ’ ) ; 

for  k3  = 1:0+1 

disp(sprintf  ( ’Row  block  7.d  of  7,d  rows  of  H3\n’  ,k3,0+l) ) ; 

Tempi  = vdp_interp(0,t0(k3) , Weightsl) ; 
for  j3  = 1 :P+1 

Weights2  = colloc(sP(j3)+omega,2*pi) ; 
if  ( j 3 ==  1) 

Temp4  = (Templ*ImZ_inv_Z  + Tempi) *B_ cell {1 , 1}; 

Temp4  = w(l) *Temp4 ; 

H3(2*k3-1 , 1)  = Temp4(l , 1) ; 

H3(2*k3-1 ,2)  = Temp4(l ,2) ; 

H3 (2*k3 , 1)  = Temp4 (2,1) ; 
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H3(2*k3,2)  = Temp4 (2 , 2) ; 

else 

if  (t0(k3)  < sP(j3)  + omega) 

Temp2  = zeros (2, 2); 
elseif  (t0(k3)  ==  sP(j3)  + omega) 

Temp2  = eye (2); 

else 

Temp2  = vdp_interp(sP(j3)+omega,tO(k3) ,Weights2) ; 
Temp3  = vdp_interp(sP(j3)+omega,2*pi ,Weights2) ; 
Temp4  = (Tempi *ImZ_inv*Temp3  + Temp2) *B_cell{j3 , 1} ; 
Temp4  = w(j3)*Temp4; 

H3(2*k3-1 ,2*j3-l)  = Temp4(l,l); 

H3(2*k3-l,2*j3)  =Temp4(l,2); 

H3 (2*k3 , 2*j  3-1)  = Temp4(2,l); 

H3(2*k3,2*j3)  = Temp4(2,2); 

end 

end 

clear  Weights2; 

end 

end 

/(Computing  H4 
disp( ’Computing  H4’); 
for  14  = 1:0+1 

disp (sprintf  ( ’ Row  block  %d  of  °/0d  rows  of  H4\n’ ,14,0+1) ) ; 
Weights3  = colloc (t0(14) ,2*pi) ; 
if  (14  ==  1) 

Tempi  = vdp_interp(0 ,2*pi .Weightsl) ; 
for  k4  = 1:0+1 

Temp2  = vdp_interp(0,t0(k4) , Weightsl) ; 

Temp3  = Temp2*ImZ_inv*Templ  + Temp2; 

Temp3  = u(l)*Temp3; 

H4(2*k4-l,l)  = Temp3(l , 1) ; 

H3 (2*k4-l , 2)  = Temp3(l ,2) ; 

H3(2*k4, 1)  = Temp3 (2 , 1) ; 

H3(2*k4,2)  = Temp3 (2 , 2) ; 

end 

else 

Tempi  = vdp_interp(t0(14) ,2*pi ,Weights3) ; 
for  k4  = 1:0+1 

if  (tO  (14)  > tO (k4) ) 

Temp2  = zeros(2,2); 
elseif  (t0(14)  ==  t0(k4)) 

Temp2  = eye (2); 

else 

Temp2  = vdp_interp(t0(14) ,t0(k4) ,Weights3) ; 

Temp3  = vdp_interp(tO (14) , 2*pi , Weights3) ; 

Temp4  = vdp_interp(0,t0(k4) .Weightsl) ; 

Temp5  = Temp4*ImZ_inv*Temp3  + Temp2; 

Temp5  = u(14)*Temp5; 

H4(2*k4-1, 2*14-1)  = Temp5(l,l); 

H4(2*k4-l,2*14)  =Temp5(l,2); 

H4(2*k4, 2*14-1)  =Temp5(2,l); 

H4 (2*k4 , 2*14)  = Temp5(2,2); 
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end 

end 

end 

clear  Weights3; 

end 

’/.Computing  final  bound 

H5  = H3*pinv(IP  - H1)*H2  + H4; 

M = norm(H5 , inf) ; 

fprintf (f id , ’ \n\nPrinting  M bound\n’); 
fprintf  (fid,  ’’/.15.8e\n’  ,M)  ; 
disp( ’Final  M bound’) 

M 

°/,disp( ’Press  Enter  to  Continue’); 

’/.pause ; 

’/.Final  error  estimates 

abs_omega  = abs(a_bar(l)) ; 

disp( ’ lambda_0 ’ ) ; 

lambda_0  = sqrt (2*pi) *abs_alpha 

fprintf (fid ,’ \n\n  lambda_0  = ’); 

fprintf  (f  id , ’ °/„15 . 8e\n  ’ , lambda_0)  ; 

disp( ’ om_min_4L0r  must  be  >=  O’); 

om_min_4L0r  = omega-4*lambda_0*norm_r 

fprintf (fid ,’ \n\n  om_min_4L0r  = ’); 

fprintf  (f  id , ’ °/»15 . 8e\n’  , om_min_4L0r)  ; 

lambda_l  = M* (1+sqrt (2*pi) *abs_alpha*norm_J) 

fprintf (fid, ’ \n\n  lambda_l  = ’); 

fprintf  (f  id , ’ °/»15 . 8e\n’  , lambda_l)  ; 

lambda_2  = (lambda_l/ (abs_omega*M) ) * ( l+2*M*cal_B) 

fprintf (fid ,’ \n\n  lambda_2  = ’); 

fprintf  (f  id , ’ °/.15 . 8e\n’  , lambda_2)  ; 

Ell  = 32*cal_K*lambda_l~2 ; 

E12  = (16*lambda_0*lambda_l/abs_omega)*(abs_omega+2*lambda_0*norm_r) ; 

E13  = 2*lambda_0'‘2  * (cal_K*norm_abs_dxt~2  + cal_B*norm_abs_ddxt)  ; 

E14  = 8*cal_B*lambda_0*lambda_2 ; 
disp(’El  must  be  <=  1’) 

El  = (Ell  + E12  + E13  + E14)*norm_r 
fprintf (fid, ’\n\n  El  = ’); 
fprintf  (fid,  ’°/.15.8e\n’  ,E1)  ; 

E21  = 32*cal_K*lambda_l~2; 

E22  = 4*lambda_0*lambda_l* (2*cal_K*norm_abs_dxt  + cal_B/abs_omega) ; 

E23  = (16*lambda_l"2*lambda_0/3)*((16/abs_omega)+4)*norm_r ; 

E24  = 4*lambda_l*lambda_0* (2*cal_K*norm_abs_dxt  + cal_B*(l  +2/abs_omega) ) ; 
E25  = 8*cal_B*lambda_2*lambda_0/abs_omega ; 

E26  = 256*cal_K*lambda_l*lambda_2*lambda_0*norm_r/ (3*abs_omega) ; 

E27  = 2*cal_B*norm_abs_ddxt*lambda_0~2 ; 
disp(’E2  must  be  strictly  < 1’) 

E2  = (E21  + E22  + E23  + E24  + E25  + E26  + E27)*norm_r 
fprintf (fid ,’ \n\n  E2  = ’); 
fprintf  (f id , ’ 7.15 . 8e\n’  ,E2)  ; 
disp( ’Final  errors’) 

Final_error_x_star  = 4*lambda_l*norm_r 
Final_error_omega_star  = 2*lambda_0*norm_r 
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fprintf (f id, ’\n\n  Final_error_x_star  = J); 
f print f (f id , ’ %15 . 8e\n’ , Final .err or _x_ star) ; 
fprintf (fid , ’ \n\n  Final_error_omega_star  = ’); 
fprintf (fid, ’%15.8e\n’ , Final_error_omega_star) ; 

f close (f id) ; 


18.1.  colloc.m. 

function  Weights  = colloc(a,b) 

global  D_hat 

global  A_hat 

global  a_bar 

global  NC 

global  zj 

global  lambda 

global  m 


Q = f ix( (b-a)/a_bar(l , 1))+1 ; 

B = zeros(2,2, 1) ; tj  = zeros (NC+1 , 1) ; g = zeros(2,l); 

WO  = zeros (2* (NC+1) , 1) ; Weights  = zeros (2* (NC+1) , 2 , Q) ; 

B_hat  = zeros (2* (NC+1) ,2*(NC+1) ,Q) ; 

%solve  for  weights  on  first  step  interval 
omega  = a_bar(l,l); 

M = ((D_hat  - (omega/2) *A_hat) “ (-1) ) ; 

WO (NC+1 , 1)  = 1;  WO (2* (NC+1) , 1)  = 0; 

Weights (:, 1 , 1)  = M*W0; 

WO (NC+1 , 1)  = 0;  WO (2* (NC+1 ) , 1)  = 1; 

Weights (: ,2,1)  = M*W0; 

"/other  intervals 
for  i = 2:Q 

B_hat (NC+1 , 1 , i)  = 2/omega;  B_hat (2* (NC+1) ,NC+2 , i)  = 2/omega; 
tj  = (omega/2) *zj  + (l/2)*(2*a  + (2*i-l) *omega) ; 
for  k = 1 : length(t j )-l 

g(l,l,l)  — a_bar (2 , 1 , 1) *cos (t j (k) ) ; 
g(2,l,l)  = -a_bar (2 , 1 , 1) *sin(t j (k) ) ; 
for  n = 2:m 

g(l,l,l)  = g(l,l,l)  + a_bar (n*2, 1) *cos (n*t j (k) ) + a_bar (n*2-l , l)*sin(n*t j (k) ) ; 
g (2 , 1 , 1)  = g(2,l,l)  - n*a_bar (n*2 , 1) *sin(n*t j (k) ) + n*a_bar (n*2-l , 1) *cos (n*t j (k)) 

end 

B(2,l,l)  = (~2*lambda/a_bar (1 , 1 , 1) )*g(l , 1 , 1) . *g(2, 1 , 1) ; 

B(2,2,l)  = (lambda/a_bar(l,l))*(l-g(l ,1,1) . ~2) ; 

"/.fill  up  the  B_hat_i  matrix 
B_hat(k,k,i)  = B(l,l,l) ; 

B_hat(k, (NC+l)+k,i)  = B(l,2,l); 

B_hat ( (NC+1) +k,k, i)  = B(2,l,l); 

B_hat ( (NC+l)+k, (NC+l)+k,i)  = B(2,2,l); 

end 

Weights ( : ,l,i)  = M* (omega/2) *B_hat (1 : 2* (NC+1) ,1:2* (NC+1) , i) ^Weights ( : ,l,i-l) ; 

Weights ( : ,2,i)  = M* (omega/2) *B_hat (1 : 2* (NC+1) ,1:2*(NC+1) , i) ^Weights ( : ,2,i-l) ; 

end 


18.2.  Derivative-series. m. 
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function  g=Derivative_series (t , a_bar ,m) 
g(l , : ) = -a_bar(2, l)*sin(t) ; 
g(2,:)  = -a_bar(2, l)*cos(t) ; 
for  n = 2:m 

g(l,:)  = g C 1 , : ) ~ n*a„bar (n*2 , 1) *sin(n*t)  + n*a_bar (n*2-l , 1) *cos (n*t) ; 
g(2,:)  = g(2,:)  - (n~2) *a_bar (n*2 , 1) *cos (n*t)  - (n"2) *a_bar (n*2-l , 1) *sin(n*t) ; 

end 

18.3.  Galerkin_series.m. 

function  g=Galerkin_series (t ,a_bar ,m) 
g = a_bar (2 , 1) *cos (t) ; 
for  n = 2:m 

g = g + (a_bar(n*2, l)*cos(n*t)  + a_bar (n*2~l , 1) *sin(n*t) ) ; 

end 

18.4.  lagrint.m. 

function  yi  = lagrint(x,y,xi) 

7.  lagrint  Interpolation  with  Lagrange  polynomials  of  arbitrary  degree 

•/. 

7.  Synopsis:  yi  = lagrint (x,y,xi) 

1 

70  Input:  x,y  = tabulated  data 

7.  xi  = point  where  interpolation  is  to  be  evaluated 

1 

7,  Output:  yi  = value  of  y at  x = xi  obtained  via  interpolation  with 

7.  polynomial  of  degree  n-1,  where  length(y)  = length(x)  = n 

dxi  = xi  - x;  7,  vector  of  xi  - x(l),  xi  - x(2)  , ...  values 

n = length(x)  ; 7o  degree  of  polynomial  is  n-1 

L = zeros  (size  (y) ) ; 7.  preallocate  L for  speed 

7,  Refer  to  section  10.2.2  in  text  for  explanation  of  vectorized  code 
7,  used  to  compute  Lagrange  basis  functions,  L(j) 

L(l)  = prod  (dxi  (2  :n) ) /prod(x(l)  -x(2  :n) ) ; 7«  j = 1 

L(n)  = prod(dxi(l  :n-l)  )/prod(x(n)-x(l  :n-l) ) ; 7»  j = n 
for  j=2:n-l 

num  = prod(dxi (1 : j -1) )*prod(dxi ( j+1 :n) ) ; 
den  = prod(x( j ) -x(l : j -1) ) *prod(x( j ) -x( j +1 :n) ) ; 

L(j)  = num/den; 
end 

yi  = sum(y.*L);  7.  Evaluate  Polynomial:  sum  of  y(j)*L(j),  j=l..n 


18.5.  Phase_series.m. 

function  x_hat=Phase_series (t , a_bar ,m) 
x_hat ( : , l)=a_bar (2 , l)*cos(t) ; 
x_hat( : ,2)=-a_bar(2, l)*sin(t) ; 
for  n = 2:m 

x_hat(:,l)  = x_hat(:,l)  + a_bar (n*2 , 1) *cos  (n*t)  + a_bar(n*2-l , l)*sin(n*t)  ; 
x_hat(:,2)  = x_hat(:,2)  - n*a_bar (n*2, l)*sin(n*t)  + n*a_bar (n*2-l , 1) *cos (n*t) ; 


end 
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18.6.  Van_der_Pol.m. 
function  vdp  = Van_der_Pol (A) 

global  m N CS  MO  M2  VO  VI  V2  T lambda  DM 

vdp  = (A (1) "2) *xdd(A)  + A( 1) *lambda* (delay_x(A) . ~2  - 1) . *delay_xd(A)  + x(A) ; 

18.7.  vdp_interp.m. 

function  Z = vdp_interp (startt , eval_pt .Weights) 

%VDP_INTERP  Lagrange  interpolation  using  Chebyshev  Points  to  generate  a 2 x 

%2  matrix  Z of  the  fundamental  solution  for  Van  der  Pol  equation.  Frequency  in  a_bar(l) . 

global  a_bar 

global  NC 

global  zj 

omega  = a_bar(l,l); 

kk  = fix((eval_pt  - startt) /omega  )+l; 

zz  = (2/omega) *eval_pt  - (2*startt  + (2*kk-l) *omega) /omega  ; 
yy(l :NC+1 ,1,1)  = Weights ( 1 : NC+1 , 1 ,kk) ; 

Z(l,l)  = lagrint (zj ,yy ,zz) ; 

yy (1 :NC+1 ,1,1)  = Weights (NC+2 : 2* (NC+1) , 1 ,kk) ; 

Z (2 , 1)  = lagrint (zj ,yy,zz) ; 

yy(l :NC+1 ,1,1)  = Weights(l :NC+1 ,2,kk) ; 

Z(l,2)  = lagrint (zj , yy , zz) ; 

yy(l :NC+1 ,1,1)  = Weights(NC+2 : 2* (NC+1) ,2 , kk) ; 

Z (2 , 2)  = lagrint (z j ,yy,zz) ; 

18.8.  Vdp_series.m. 

function  x_hat=Vdp_series (t , a_bar ,m) 

%set  up  for  dde23 
x_hat (1 , : )=a_bar (2 , l)*cos(t) ; 
x_hat (2 , : )=-a_bar (2 , l)*sin(t) ; 
for  n = 2:m 

x_hat(l,:)  = x_hat(l,:)  + a_bar (n*2 , 1) *cos (n*t)  + a_bar (n*2-l , 1) *sin(n*t) ; 
x_hat(2,:)  = x_hat(2,:)  - n*a_bar (n*2 , 1) *sin(n*t)  + n*a_bar (n*2-l , 1) *cos (n*t) ; 

end 
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